Gradient-based approximate design optimization

Pełen tekst

(1)Gradient-based Approximate Design Optimization.

(2)

(3) Gradient-based Approximate Design Optimization. Proefschrift. ter verkrijging van de graad van doctor aan de Technische Universiteit Delft, op gezag van de Rector Magnificus prof.dr.ir. J.T. Fokkema, voorzitter van het College voor Promoties, in het openbaar te verdedigen. op dinsdag 1 november 2005 te 13.00 uur. door. Koen VERVENNE. werktuigbouwkundig ingenieur geboren te Wemeldinge.

(4) Dit proefschrift is goedgekeurd door de promotor: Prof.dr.ir. A. van Keulen. Samenstelling promotiecommissie: Rector Magnificus, Prof.dr.ir. A. van Keulen, Prof.dr. V.V. Toropov, Prof.dr. Z. Gürdal, Prof.dr. R.T. Haftka, Prof.dr.ir. D. den Hertog, Prof.dr.ir. A. de Boer, Dr.ir. L.F.P. Etman,. voorzitter Technische Universiteit Delft, promotor University of Bradford Technische Universiteit Delft University of Florida Universiteit van Tilburg Universiteit Twente Technische Universiteit Eindhoven. Published and distributed by: DUP Science DUP Science is an imprint of Delft University Press P.O. Box 98 2600 MG Delft The Netherlands Telephone: +31 15 27 85 678 Telefax: + 31 15 27 85 706 E-mail: info@library.tudelft.nl ISBN 90-407-2608-6 Keywords: Gradient, Response surface, Design optimization c Copyright 2005 by K. Vervenne All rights reserved. No part of the material protected by this copyright notice may be reproduced or utilized in any form or by any means, electronic or mechanical, including photocopying, recording or by any information storage and retrieval system, without written permission from the author..

(5) Contents 1. Introduction 1.1 Link to the ADOPT project . . . . . . . . . . . . . . . . . . . . . 1.2 Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 2. Accuracy Improvement of Semi- Analytical Design Sensitivities by Laplacian Smoothing 5 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.2 Semi-analytical design sensitivities . . . . . . . . . . . . . . . . . 7 2.3 Refined semi-analytical design sensitivities . . . . . . . . . . . . 8 2.4 Laplacian smoothing . . . . . . . . . . . . . . . . . . . . . . . . 9 2.5 Laplacian smoothing with remainder . . . . . . . . . . . . . . . . 10 2.6 Analytical beam example . . . . . . . . . . . . . . . . . . . . . . 12 2.7 Numerical examples . . . . . . . . . . . . . . . . . . . . . . . . 15 2.7.1 Rectangular strip . . . . . . . . . . . . . . . . . . . . . . 15 2.7.2 Cylindrical panel . . . . . . . . . . . . . . . . . . . . . . 21 2.8 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24. 3. An Alternative Approach to Response Surface Building using Gradient Information 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Response surface building . . . . . . . . . . . . . . . . . . . . . 3.3 Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.1 Example I . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.2 Example II . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 27 27 29 33 33 40 43. Gradient-Enhanced Response Surface Building 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 GERS building techniques . . . . . . . . . . . . . . . . . . . . . 4.2.1 Weighted Least Squares . . . . . . . . . . . . . . . . . .. 45 45 47 48. 4. v. 1 3 3 4.

(6) vi. 4.3. 4.4 5. 4.2.2 Two-step blending . 4.2.3 Pareto-based GERS . Case Studies . . . . . . . . . 4.3.1 Example 1 . . . . . 4.3.2 Example 2 . . . . . 4.3.3 Example 3 . . . . . Discussion and Conclusions. . . . . . . .. . . . . . . .. . . . . . . .. . . . . . . .. . . . . . . .. . . . . . . .. . . . . . . .. Gradient-Enhanced Radial Basis Functions 5.1 Introduction . . . . . . . . . . . . . . . . 5.2 Functions with different length-scales . . 5.3 Radial Basis Functions . . . . . . . . . . 5.4 Gradient enhanced radial basis functions . 5.5 Examples . . . . . . . . . . . . . . . . . 5.5.1 1-dimensional Griewank function 5.5.2 N-dimensional Griewank function 5.6 Discussion and Conclusions . . . . . . . 5.7 Acknowledgement . . . . . . . . . . . .. . . . . . . .. . . . . . . . . .. . . . . . . .. . . . . . . . . .. . . . . . . .. . . . . . . . . .. . . . . . . .. . . . . . . . . .. . . . . . . .. . . . . . . . . .. . . . . . . .. . . . . . . . . .. . . . . . . .. . . . . . . . . .. . . . . . . .. . . . . . . . . .. . . . . . . .. . . . . . . . . .. . . . . . . .. . . . . . . . . .. . . . . . . .. . . . . . . . . .. . . . . . . .. . . . . . . . . .. . . . . . . .. 49 50 53 53 61 69 76. . . . . . . . . .. 83 84 86 87 89 90 90 91 92 94. 6. Efficiency Improvement of Response Surface Building using Fast Reanalysis Methods 95 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 6.2 Response surface method . . . . . . . . . . . . . . . . . . . . . . 97 6.3 Fast reanalysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 6.3.1 Implementation details . . . . . . . . . . . . . . . . . . . 99 6.3.2 Accuracy problems for shape design variables . . . . . . . 100 6.3.3 Region of accurate results . . . . . . . . . . . . . . . . . 101 6.4 Response surface building using fast reanalysis . . . . . . . . . . 101 6.4.1 Design of experiments . . . . . . . . . . . . . . . . . . . 103 6.4.2 Error analysis . . . . . . . . . . . . . . . . . . . . . . . . 104 6.5 Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 6.5.1 Truss example . . . . . . . . . . . . . . . . . . . . . . . 105 6.5.2 Strip example . . . . . . . . . . . . . . . . . . . . . . . . 106 6.6 Discussion and Conclusions . . . . . . . . . . . . . . . . . . . . 111. 7. Sequential Approximate Optimization of Composite Laminates 7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2 Laminate failure analysis . . . . . . . . . . . . . . . . . . . 7.3 Optimization problem formulation . . . . . . . . . . . . . . 7.4 Optimization approach . . . . . . . . . . . . . . . . . . . . 7.5 Example I . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . .. . . . . .. . . . . .. 113 113 114 117 117 118.

(7) vii 7.6 7.7 8. Example II . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 Discussion and Conclusions . . . . . . . . . . . . . . . . . . . . 120. Discussion and future work 125 8.1 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 8.2 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126.

(8) viii.

(9) Chapter 1 Introduction In order to improve product and process design, optimization techniques are frequently used. The initial academic interest in the topic is now also present in industry. For example, in the automotive and aerospace industry, optimization has become a common design tool. The optimization process usually consists of finding a design with better features as compared to the original design, but which still fulfills the requirements. In order to carry out an optimization task, variables that affect the performance of the design have to be selected. These variables are called design variables. It is also necessary to exactly define what has to be optimized. A single or multiple objective functions have to be formulated in terms of the design variables. Finally, the requirements put on the performance of the design have to be defined. These requirements are called constraint functions and they also depend, generally implicitly, on the design variables. For example, in case of structural optimization, the dimensions of a structure and the material properties might be selected as design variables. The weight of the structure might be chosen as the objective function. Examples of constraints are maximum deflection or maximum equivalent stress. The evaluation of the objective functions and constraint functions for a particular design can be very time consuming. For example, in structural optimization the finite element method is often used, which leads to large computation times. More often, multiple disciplines are involved in calculation of the objectives and constraints. This may lead to an even larger computation time. By means of parallel computing on large clusters of computers this computation time can be reduced. Despite developments in computer technology during the last decades, the demand for computational power is often still larger than the available resources. Another possibility to reduce the computational load of an optimization process is to approximate the numerically expensive objective and constraint functions. A popular technique to construct these approximations is the response sur1.

(10) 2. Chapter 1. face method (Myers and Montgomery 1995; Khuri and Cornell 1996). Based on a limited number of expensive function evaluations, approximation functions are created that can be used during the entire optimization process or during one or more cycles of the optimization process. The response surface method is considered to be a method to reduce the so-called curse of dimensionality, which refers to the enormous increase of the number of function evaluations if the number of design variables is increased. In some applications, in addition to the usual function evaluations, sensitivity information can be calculated inexpensively. The sensitivities provide the derivatives of the response functions with respect to the design variables. In (van Keulen, Haftka, and Kim 2005) a review of options for structural sensitivity calculations can be found. The construction of response surface approximations might benefit from sensitivity information. If sensitivities are used, the number of expensive function evaluations could be reduced while the accuracy of the approximation remains sufficient. Another possibility is that the number of design variables could be increased by using sensitivity information. In this way, the sensitivity information can help to solve problems with a larger number of design variables and thus reduce ”the curse of dimensionality”. In literature several methods have been proposed in order to include sensitivity information into the response surface. However, these methods are not well crystallized and it is not clear which approach is to be preferred. In the present thesis several methods for gradient-enhanced response surface building have been studied. Simple numerical examples have been used in order to compare the various methods. In case of structural finite element analysis, the semi-analytical method is often used to compute the sensitivities. Calculation of semi-analytical sensitivities only requires a fraction of the time of the original analysis. However, in case of shape design variables the method leads to severe accuracy problems (Barthelemy, Chon, and Haftka 1988; Pedersen, Cheng, and Rasmussen 1989; de Boer and van Keulen 1997). In the present thesis it has been investigated whether these accuracy problems can be reduced by using Laplacian smoothing. Instead of using sensitivity information to build the response surface, additional function evaluations, calculated by means of fast reanalysis, could be used. A fast reanalysis technique has been developed by Kirsch (Kirsch 2002; Kirsch 2003) and provides a way to cheaply calculate the response of a structure for a modified design. By using fast reanalysis in combination with the response surface method, the efficiency of the response surface building could be improved. Simple finite element problems have been used to demonstrate this approach. Finally, as a more concrete problem from practice, the present thesis will focuss on optimization of composite laminates. This problem has been formulated by the company Airborne Development and contains several complications that prevent straightforward solution by means of standard optimization techniques..

(11) Introduction. 3. In the present problem uncertainties in the design variables play a role, as well as discrete design variables and discontinuities. In this thesis, the emphasis is on the discontinuities in the response function, and it will be investigated whether a sequential approximate optimization method can be used to solve the optimization problem.. 1.1. Link to the ADOPT project. The present research has been carried out as part of the ADOPT project (Approximate Design OPTimization). This project is a collaboration between Delft University of Technology and Eindhoven University of Technology. The objective of the ADOPT project is to develop, implement and test computer-aided optimization strategies based on the sequential approximate design concept. These optimization strategies should be able to deal with combinations of discontinuities, uncertainties, and discrete design variables. Applications used in the ADOPT project are optimization of composite structures, MEMS (Micro Electrical Mechanical Systems), dynamical systems and manufacturing systems. The application used in the present project focuses on optimization of composite structures.. 1.2. Objectives. The research presented in this thesis has the following objectives: • Development of accurate sensitivity analysis based on the semi-analytical method. This method is known for severe inaccuracies for shape design variables. The influence of the mesh perturbation method on the accuracy of the sensitivities will be investigated. • Development of response surface techniques based on both function values and gradients. Several methods to include gradient information into the response surface will be formulated. Their derivation will be presented and advantages and disadvantages will be discussed based on several examples. • Application of the developed techniques to both analytical and numerical examples. The optimization of composite tubes will be used as a representative example from practice..

(12) 4. 1.3. Chapter 1. Outline. The outline of the thesis is as follows. Chapter 2 discusses the use of Laplacian smoothing in order to improve the accuracy of semi-analytical sensitivities. Chapter 3 explains how these sensitivities can be used to construct response surfaces. Here, a two-step approach to construct the response surface will be presented. Construction of gradient-enhanced response surfaces can be considered as a multi-objective problem which will be shown in Chapter 4. The use of radial basis functions in order to construct response surfaces based on function values and gradients will be presented in Chapter 5. Chapter 6 focuses on efficiency improvement of response surface building using fast reanalysis methods. In Chapter 7 the optimization of composite laminates using a sequential approximate optimization approach will be discussed. Finally, Chapter 8 contains a discussion and general conclusions. Chapters 2 to 6 are reproductions from papers presented at conferences and journal publications. References to these papers are given at the first page of the corresponding chapter. As a consequence of the fact that Chapters 2 to 6 are reproductions, some information in this thesis is superfluous. On the other hand, it enables selective reading of each of the chapters..

(13) Reproduced from: K. Vervenne and F. van Keulen (2001), Accuracy improvement of semi-analytical design sensitivities by Laplacian smoothing, 42nd Structures, Structural Dynamics, and Materials Conference and Exhibit, 16-19 April, Seattle, WA, USA, AIAA-2001-1497. Chapter 2 Accuracy Improvement of SemiAnalytical Design Sensitivities by Laplacian Smoothing Abstract Sensitivity analysis by means of the semi-analytical method often results in severe accuracy problems for shape design variables. It turns out that the method used to generate the perturbed meshes has an important effect on the accuracy of the design sensitivities. In the present paper it will be shown that the accuracy of the design sensitivities can be improved by using Laplacian smoothing instead of a boundary node approach. When Laplacian smoothing is used design perturbations are distributed more equally over the finite element mesh. Implementation of the method is easy in an existing preprocessor. The proposed method is applied to semi-analytical and refined semi-analytical design sensitivities. In both cases more accurate design sensitivities are obtained. The accuracy of the design sensitivities has been studied by means of an analytical beam example and several numerical examples. Key words: design sensitivities, semi-analytical, refined semi-analytical, mesh perturbation, accuracy, Laplacian smoothing.. 2.1. Introduction. The semi-analytical method is often used for computation of design sensitivities in a finite element context. This method is easy to implement in an existing finite 5.

(14) 6. Chapter 2. element code and computationally efficient. However, severe accuracy problems can occur, especially for shape design variables. Several methods to improve the accuracy of the semi-analytical method can be found in the literature. The use of higher order finite difference schemes has been proposed by Barthelemy, Chon, and Haftka (1988), Cheng, Gu, and Zhou (1989) and Olhoff and Rasmussen (1991). Olhoff and Rasmussen (1991) developed the “exact” formulation to obtain accurate sensitivities. This “exact” method exploits the specific structure of the finite element matrices. Alternative forward/backward finite difference schemes have been proposed by Cheng, Gu, and Wang (1991). Mlejnek (1992) used consistency equations for rigid body modes to increase the accuracy. Oral (1996) presents an improved semi-analytical method using von Neumann series. So-called refined semi-analytical sensitivities have been developed by van Keulen and de Boer (1998). The method was also adopted by Parente and Vaz (2001). The refined semi-analytical method is based on analytical differentiation of rigid body modes particularly the rigid body rotations. In the semi-analytical method derivatives of the stiffness matrix with respect to the design variables are approximated by finite difference schemes. This means that perturbed meshes corresponding to slightly changed values of the design variables are necessary. The generation of these meshes is given only little attention in the references mentioned above and is therefore the focus of the present study. In the present paper it will be shown that in addition to the improvements advocated in the references given above, the accuracy of the semi-analytical method can also be improved by increasing the quality of the perturbed meshes. In literature several methods can be found in order to create finite element meshes corresponding to relatively large changes of the design variables. A typical application is mesh updating during shape optimization. A method is to relocate boundary nodes only. This procedure is called the boundary node approach (Lindby and Santos 1997). The boundary node approach can result in relatively large perturbations of the boundary elements. Several methods for computation of the location of the inner nodes have been described. The so-called geometrical approach is described by Schramm and Pilkey (1993), Lindby and Santos (1997) and Hwang, Choi, and Chang (1997) and is based on an explicit relation between the nodal coordinates and the shape design variables. In the physical approach (Beckers 1991) the displacements of the boundary nodes are used as prescribed displacements in a finite element computation. Perturbed meshes can also be generated by using a spring analogy (Bugeda and Onate 1993) which is very similar to the physical approach. A disadvantage connected to both the physical and the spring model approach is that a large set of equations has to be assembled and stored subsequently. Finally, perturbed meshes can be generated by using Laplacian smoothing (Bugeda and Onate 1993). In contrast to the physical approach and the spring model, the Laplacian smoothing.

(15) Accuracy improvement of semi-analytical design sensitivities. 7. benefits from an iterative setup. Generally, the mesh perturbation techniques as sketched above are being used for the generation of meshes with relatively large design perturbations. That is, the design perturbations are the result of design adaptations being controlled, for example, by an optimization algorithm. In the present paper, however, the boundary node approach and Laplacian smoothing have been applied to the computation of semi-analytical design sensitivities. Thus, the Laplacian smoothing is used to generate meshes corresponding to relatively small design perturbations. By means of an analytical beam example and several numerical examples it will be shown that perturbation of all nodes of the mesh by means of Laplacian smoothing leads to more accurate design sensitivities as compared to a boundary node approach. The outline of the paper is as follows. First, for self-containment, the semianalytical and the refined semi-analytical method for computation of design sensitivities will be briefly discussed. Then, Laplacian smoothing and its application to design sensitivity analysis will be explained. An analytical beam example will be used in order to show that the accuracy of design sensitivities is dependent on the method used to create perturbed meshes. Thereafter, numerical examples will be presented and finally the results will be discussed.. 2.2. Semi-analytical design sensitivities. The finite element equations describing linear elastic structural behaviour are Ku = f.. (2.1). In this equation, K is the stiffness matrix, u is the vector of nodal degrees of freedom and f is the vector of nodal loads. First order sensitivities of the displacement vector are obtained by differentiation of (2.1) with respect to design variable si , which results in Ku,i = f,i − K,i u. (2.2) Here ...,i denotes differentiation with respect to si . The right-hand side of (2.2) is normally referred to as the pseudo-load vector. In the semi-analytical method derivatives of the stiffness matrix with respect to the design variables are approximated by means of finite difference schemes. In case a forward finite difference scheme is used this leads to K,i ≈. K (si + ∆si ) − K (si ) . ∆si. (2.3). The displacement sensitivities can be solved from (2.2), using the stiffness matrix decomposition which is already available from the static analysis based on (2.1)..

(16) 8. 2.3. Chapter 2. Refined semi-analytical design sensitivities. A rigorous improvement of the semi-analytical method has been developed by van Keulen and de Boer (1998) and can be used in order to improve the accuracy of the semi-analytical design sensitivities. In the so-called refined semi-analytical method the element displacement vector is decomposed into a part leading to deformations and a part leading to rigid body motions. This yields ue = uεe + αk rke. with. αk =. ue ·rke . rke ·rke. (2.4). Here subscripts e denote quantities at element level, ue is the element displacement vector, uεe refers to the deformation part of the element displacement vector and rke denotes the set of rigid body modes of the element. Note, use has been made of summation convention for repeated indices. The rigid body modes should fulfill Ke rke = 0. (2.5). and by differentiation of this equation with respect to design variable si it follows Ke ,i rke + Ke rke ,i = 0.. (2.6). These equations will be referred to as consistency equations. By means of the decomposition of the displacement vector (2.4) and the consistency equations for rigid body modes (2.6) the pseudo-load vector of a single element can be written as qe = fe ,i − Ke ,i uεe + αk Ke rke ,i .. (2.7). From this equation it can be seen that the accuracy of the pseudo-load vector is improved compared to the semi-analytical method because the derivative of the element stiffness matrix is partly replaced by analytical derivatives of the rigid body modes. In contrast to the derivatives of the element matrices (Ke ), the derivatives of the rigid body modes are easy to compute analytically. A second improvement of the pseudo-load vector is possible by a decomposition of the pseudo-load vector into a self-equilibrating part and a non-self-equilibrating part. Possible errors in the self-equilibrating part of the pseudo-load vector do not have a large effect on the accuracy of the design sensitivities. This can be concluded from Saint Venant’s principle: the effect of inaccuracies in self-equilibrating terms tends to damp out. The accuracy of the non-self-equilibrating part can be improved by using again the consistency equations and analytical differentiation of rigid body.

(17) Accuracy improvement of semi-analytical design sensitivities modes. The pseudo-load vector of a single element then becomes " # k ·K uε r e ,i e e qe = fe ,i − Ke ,i uεe + rke + rke ·rke # " rke ,i ·Ke uεe k re + αk Ke rke ,i . k k re ·re. 9. (2.8). By means of the described refined formulation a significant accuracy improvement of the semi-analytical method can by obtained. Another advantage is that the implementation in an existing finite element code is easy. Since the pseudo-load vector is more complicated compared to the standard formulation, there will be an increase in computation time which is proportional to the number of finite elements in the mesh.. 2.4. Laplacian smoothing. In this section Laplacian smoothing will be explained in more detail and its application to design sensitivity analysis will be discussed. Throughout the discussion it will be assumed that a finite element mesh consisting of triangular elements is available. Laplacian smoothing is a smoothing algorithm which modifies the positions of the internal nodes of the mesh while keeping the topology of the mesh unchanged (Vollmer, Mencl, and Muller 1999). During Laplacian smoothing the coordinates of each interior node of the mesh will be replaced by the average of the coordinates of adjacent nodes. A node is considered adjacent if the current node and the adjacent node belong at least to the set of nodes of a single element. In other words the nodes must be connected by at least one element. The position vector after each iteration is given by m. i rj ∑ j=1 ri := , mi. (2.9). where r j are the position vectors of the mi nodes connected to the i-th node. There are two variants of Laplacian smoothing. In the first variant the nodal coordinates are updated only after the averages for all internal nodes have been computed. This variant is called the simultaneous version. In the second variant, the new nodal coordinates are updated immediately. This variant is called the sequential version. A consequence of updating the nodal coordinates simultaneously is that more memory is required for storing the old positions until all new ones are computed. An advantage is that better results are obtained in the sense that the element condition is improved. However both the simultaneous and.

(18) 10. Chapter 2. the sequential version converge towards the same solution (Vollmer, Mencl, and Muller 1999). In this work the sequential version has been used. This choice was prompted because implementation of the sequential version is less involved than implementation of the simultaneous version. Laplacian smoothing is carried out using only the adjacent nodes. The central point could also be included, but this modification does not result in a better quality of the mesh. It only delays the smoothing process compared with the original method (Vollmer, Mencl, and Muller 1999). An advantage of Laplacian smoothing is that it can also be used for curved structures. In the iteration scheme, the three-dimensional coordinates then have to be replaced by the corresponding surface coordinates. At the end of the smoothing process, the surface coordinates are transformed back to three-dimensional coordinates. A disadvantage is that Laplacian smoothing is an iterative method, but usually convergence is obtained in a limited number of iterations. During the Laplacian smoothing the nodes of the mesh are handled in sequence of increasing node number. This causes the results of the smoothing to be dependent on the node numbering. In order to reduce this dependency the nodes could be handled alternately, i.e. in increasing and decreasing order. Laplacian smoothing can be used to improve the accuracy of sensitivity analysis as follows. A finite element mesh is generated by a preprocessor. Then perturbed meshes corresponding to each value of the design variables are created by means of the boundary node approach. Laplacian smoothing is used in order to increase the quality of the perturbed meshes. Derivatives of the stiffness matrix with respect to the design variables which are required for the computation of design sensitivities are computed by means of finite difference schemes based on the original and perturbed meshes. It has been shown by Lindby and Santos (1997) that the distribution of the design perturbations influences the accuracy of the design sensitivities. Since the perturbations are distributed more uniformly if Laplacian smoothing is used in addition to a boundary node approach, more accurate design sensitivities can be expected.. 2.5. Laplacian smoothing with remainder. If Laplacian smoothing is applied to the calculation of design sensitivities as described above, the following problem could occur. Application of Laplacian smoothing to the perturbed mesh can result in node perturbations that are much larger than the imposed design perturbations, which are usually very small. This will result in inaccurate finite difference approximations and hence in inaccurate design sensitivities. The easiest way to understand this problem is to observe what happens if the.

(19) Accuracy improvement of semi-analytical design sensitivities. 11. perturbation introduced becomes infinitely small. Application of the Laplacian smoothing now basically attempts to improve the original mesh, i.e. the mesh without mesh perturbation. A way to eliminate this complication is to run the Laplacian smoothing for the unperturbed mesh until full convergence is obtained. This has two main drawbacks. First this is expensive as the iteration process may call for many iterations for full convergence. Second, if mesh refinements are included in the mesh they may get affected by the smoothing algorithm. To overcome these complications, an alternative method is proposed here. In Haftka (1993) a method is presented in order to obtain accurate global finite difference sensitivities of non-linear finite element responses. The problem studied by Haftka (1993) is related to geometrically nonlinear solutions which are normally only at our disposal with a finite accuracy. That means the remainder of the governing equations is not iterated to zero entirely. To compensate for this remainder, the remainder has to be taken into account in the global finite difference scheme. The problem described above for Laplacian smoothing can be solved by using a similar approach. In order to obtain accurate sensitivities with only a limited number of iterations, the original Laplacian smoothing algorithm can be modified as follows. First, for the original mesh the modifications that would result from one Laplacian smoothing step are computed and stored. These modifications will be called the “remainders” Ri and are defined by Ri = r0i −. 1 mi. mi. ∑ r0j ,. (2.10). j=1. where r0j are the position vectors of the mi nodes adjacent to the i-th node of the original mesh. During Laplacian smoothing of the perturbed mesh, each position vector is corrected by adding the corresponding remainder, hence 1 ri := mi. mi. ∑ r j + Ri .. (2.11). j=1. By using this approach large changes in the perturbed mesh compared to the size of the design perturbation can be prevented. This will result in more accurate finite difference approximations and consequently in more accurate design sensitivities. In order to get a better understanding of the proposed method, it is illustrative to study the situation for which the design perturbations are all kept zero. In that case r j in the righthand side of (2.11) should be replaced by r0j . Substitution of (2.10) into (2.11) now clearly shows that the nodal locations are not affected by the iteration scheme..

(20) 12. Chapter 2. There is also another interpretation possible. For that purpose we introduce the nodal perturbation, which is for node i defined as pi = ri − r0i .. (2.12). If we combine (2.10) and (2.11) and substitute (2.12), then it follows pi =. 1 mi. mi. ∑ p j.. (2.13). j=1. Hence the method of introducing the remainder is entirely equivalent to application of Laplacian smoothing to the nodal perturbations. One may question whether there is a need for two formulations, i.e. one based on the introduction of a remainder, i.e. using (2.10) and (2.11) and one based on Laplacian smoothing of the perturbations, i.e. (2.12) and (2.13). It turns out that the formulations, though equivalent, may be different from the point of implementation. Thus, depending on the software implementation at hand it may be advantageous to introduce a remainder or to work with perturbations solely. Another advantage of using a remainder or Laplacian smoothing of the perturbations is related to the presence of local mesh refinements in the mesh. If the original formulation of Laplacian smoothing would be used these local mesh refinements would be smoothed out and disappear completely. By using the alternative schemes the mesh refinements remain unaffected and only the mesh perturbation is smoothed over the mesh.. 2.6. Analytical beam example. The effect of the mesh perturbation method on the accuracy of semi-analytical design sensitivities can be demonstrated by means of an analytical beam example (Figure 2.1). This beam example has been studied previously by Barthelemy, Chon, and Haftka (1988), Pedersen, Cheng, and Rasmussen (1989), Olhoff and Rasmussen (1991), and de Boer and van Keulen (1997). For a standard semianalytical formulation the relative error in the first order design sensitivity of the tip displacement with respect to the length of the beam is dependent on the number of elements squared (de Boer and van Keulen 1997): 5 1 1 2 + 5η + 2η2 2 ε=− η η n − . 2 (1 + η)3 2 (1 + η)3. (2.14). Here η denotes the global relative design perturbation (∆L/L) and n the number of elements. It should be noted that in this case all elements of the beam have.

(21) Accuracy improvement of semi-analytical design sensitivities. 13. Figure 2.1: Beam consisting of n uniformly distributed elements and loaded by a moment at the tip.. been perturbed uniformly in order to create the perturbed mesh. This perturbed mesh could have been created by means of Laplacian smoothing. For reasons of efficiency, perturbed meshes are often created by perturbing boundary nodes only (Lindby and Santos 1997). For the beam example this corresponds to a perturbation of the tip element only. The error in the first order tip displacement sensitivity can be derived as follows. First order sensitivities of the displacement vector can be obtained from (2.2). In this case the derivative of the load vector is zero (f,i = 0). The derivative of the stiffness matrix is approximated by a finite difference scheme. Only the last element of the beam has a contribution to the error in the pseudo-load vector because only the last element is perturbed. This error can be given by eqe = −nEef ue .. (2.15) f. In this equation ue is the element displacement vector of the tip element and Ee is the error in the derivative of the stiffness matrix of an element. Using beam theory, the displacement vector for the tip element can be calculated:. ue T =. Ml l(n − 1)2 2l0 (n − 1) ln2 2l0 n . 2EI. (2.16). Here l is the element length and l0 is an arbitrary reference length. The error in the derivative of the element stiffness matrix is . 36β1  EI 12ξβ 2 Eef = − 4  nl  −36β1 12ξβ2. 12ξβ2 −36β1 4ξ2 β3 −12ξβ2 −12ξβ2 36β1 2ξ2 β3 −12ξβ2.  12ξβ2 2ξ2 β3  . −12ξβ2  4ξ2 β3. (2.17).

(22) 14. Chapter 2. In this equation ξ is defined by l/l0 and the coefficients βi are defined by 6 + 8η∗ + 3η∗2 ; 3(1 + η∗ )3 3 + 2η∗ β2 = −η∗ ; 2(1 + η∗ )2 1 , β3 = −η∗ (1 + η∗ ) β1 = −η∗. (2.18). which are dependent on the local relative design perturbation, η∗ = ∆L/l. Substitution of (2.16) and (2.17) in (2.15) gives for the error in the pseudo-load vector   3n(−6β1 n + 3β1 + 4β2 n − 2β2 ) nM  −L/l0 (6β2 n − 3β2 − 3β3 n + 2β3 )  . eqe = 2 2  (2.19) L  −3n(−6β1 n + 3β1 + 4β2 n − 2β2 )  −L/l0 (6β2 n − 3β2 − 3β3 n + β3 ) From the error in the pseudo-load vector the relative error in the first order derivative of the tip displacement can be computed by using beam theory which results in ε=−. 1 η 2 (1 + ηn)3 2 η2 + 6 n3 + 6 (η − 1) n2 − ηn + 1 .. (2.20). It should be noted that in this equation η denotes the global relative design perturbation (∆L/L) which is not equal to the local relative design perturbation (∆L/l). For the case that the beam consists of a single element (n = 1), equation (2.14) and equation (2.20) become equal, which can be checked easily. It can be seen from equation (2.20) that the error in the first order design sensitivity is dependent on n3 in case a boundary node approach is used. If a uniform perturbation is used the error is dependent on n2 (equation (2.14)). Another interesting comparison between equation (2.14) and equation (2.20) can be made if the global perturbation is replaced by the local perturbation. In case of a uniform perturbation the global perturbation is equal to the local perturbation of all elements and equation (2.14) remains the same. However, if only the tip element is perturbed the global perturbation is related to the local perturbation as follows ∆L ∆L 1 ∗ η= = = η , (2.21) L nl n where η∗ denotes the local perturbation. If this relation is substituted into equation (2.20) the error in the sensitivity becomes dependent on n2 . From these results it.

(23) Accuracy improvement of semi-analytical design sensitivities. 15. Figure 2.2: Mesh of the strip, boundary conditions and loading (150 elements).. follows that the error is dependent on n2 for both a boundary node approach and a uniform perturbation, provided the same local perturbation is used. By means of this simple example it has been shown that the mesh perturbation method has a large effect on the accuracy of the design sensitivities, provided that the global perturbation is equal. Perturbation of all elements of the mesh (Laplacian smoothing) results in more accurate sensitivities as compared to a perturbation of the last element only (boundary node approach).. 2.7. Numerical examples. The effect of Laplacian smoothing on the accuracy of the semi-analytical and refined semi-analytical design sensitivities will now be studied by means of several numerical examples.. 2.7.1 Rectangular strip As a numerical example a rectangular strip is studied. The strip is clamped at one side and loaded by a distributed bending moment at the other side. The mesh, the boundary conditions and the loading are presented in Figure 2.2. The strip is modeled using triangular shell elements (van Keulen and Booij 1996) with a thickness of 1 mm. Young’s modulus is 2 × 105 N/mm2 and Poisson’s ratio is 0.3. The design variable is chosen as the length L of the strip and its value is 100 mm. It is important to note, that the finite element model will give exact results for this problem. Thus a mesh refinement or a mesh perturbation will never affect the accuracy of the finite element model. This numerical example is used to investigate the effect of Laplacian smoothing on the accuracy of the design sensitivities. First order design sensitivities of the tip displacement with respect to the length of the strip have been computed by means of the semi-analytical and the refined semi-analytical method. The analytical solution for the tip displacement can be obtained easily from plate theory and is equal to ML2 1 − ν2 w=6 . (2.22) Eh3.

(24) 16. Chapter 2. Figure 2.3: The perturbed mesh with positive perturbation (150 elements).. Figure 2.4: The perturbed mesh with positive perturbation, after Laplacian smoothing (150 elements).. The first order design sensitivity is obtained by differentiation of the displacement with respect to the length L of the strip dw 2 = w. dL L. (2.23). In the sequel, logarithmic design sensitivities will be used, defined by w0l =. L dw = 2. w dL. (2.24). Hence the analytical logarithmic design sensitivity is equal to 2 and this value will be compared with the numerical results. Perturbed meshes have been generated by means of both the boundary node approach and Laplacian smoothing. In order to demonstrate the effect of Laplacian smoothing visually, a perturbed mesh with an extremely large mesh perturbation (η = 0.1) using the boundary node approach is depicted in Figure 2.3. The quality of the perturbed mesh has been improved by the sequential version of Laplacian smoothing. During the Laplacian smoothing process 500 iterations have been used. The perturbed mesh, after Laplacian smoothing, is shown in Figure 2.4. The convergence behavior of Laplacian smoothing applied to a mesh with a relative perturbation of 0.1 is presented in Figure 2.5. Here the maximum relative iterative perturbation, as observed during the Laplacian smoothing, is defined as ∆x = ∆X/L where ∆X is the maximum perturbation at a certain stage of the smoothing process. The largest node perturbation in y-direction is zero from the first iteration onwards. This is caused by the fact that a remainder has been used to prevent perturbations due to the initial mesh. It is seen that the perturbation converges rapidly. First order semi-analytical and refined semi-analytical design sensitivities have been computed corresponding to different values of the relative perturbation and a mesh consisting of 150 elements (Figure 2.4 and Figure 2.6). The perturbed.

(25) Accuracy improvement of semi-analytical design sensitivities 0.1. 17. x-coordinate y-coordinate. 0.08. ∆x. 0.06. 0.04. 0.02. 0 0. 2. 4. 6. 8. 10. number of iterations. Figure 2.5: Convergence of Laplacian smoothing in x- and y-direction for the perturbed mesh (mesh perturbation 0.1). The largest perturbation has been normalized with the length of the strip.. meshes have been created by using both the boundary node approach and Laplacian smoothing. During the smoothing 500 iterations have been used and the formulation with the remainder has been used. All sensitivities become inaccurate for large values of the relative perturbation because of truncation errors. For small values of the relative perturbation the results also become inaccurate because of round-off errors. If Laplacian smoothing is used, the range of relative perturbations resulting in accurate sensitivities is extended for both the semi-analytical and the refined semi-analytical method. So far only a single mesh density and a fixed number of iterations for the Laplacian smoothing has been used. This number has been selected such that convergence was guaranteed. Now the effect of the number of iterations will be studied more carefully. This will be done in conjunction with the applied design perturbation and the mesh density. For this study the following definitions will be used: the relative perturbation η = ∆L/L where ∆L is the absolute perturbation and L is the design variable, the relative mesh size h = H/L where H is the mesh size. Note that the mesh density is uniform among the mesh. The maximum relative iterative perturbation, as observed during the Laplacian smoothing, is defined as ∆x = ∆X/L where ∆X is the maximum perturbation at a certain stage of the smoothing process. For this problem the relative error in the first order refined semi-analytical design sensitivities has been computed corresponding to different numbers of iterations during Laplacian smoothing of the perturbed mesh. This relative error is defined by ε = (w0l − 2)/2,. (2.25).

(26) Chapter 2. logarithmic 1st order design sensitivity. 18. 2.06. 2.04. 2.02. 2. 1.98 -11. SA RSA SA-Laplace RSA-Laplace -10. -9. -8. -7. -6. -5. -4. -3. -2. log relative perturbation. Figure 2.6: First order design sensitivities computed by means of the semianalytical method and the refined semi-analytical method, with and without Laplacian smoothing.. where w0l is the logarithmic design sensitivity as defined before. The error has been computed for initial relative design perturbations of 10−2 and 10−4 . For these perturbations one may expect that errors in sensitivities are dominated by truncation errors rather than round-off errors. Both a coarse mesh (h = 0.04) and a fine mesh (h = 0.004) have been used. The results for the refined semi-analytical formulation are presented in Figure 2.7. The logarithm of the error has been plotted along the vertical axis. The logarithm of the maximum relative perturbation that occurred during Laplacian smoothing ∆x/η has been plotted along the horizontal axis. In all cases the error decreases if a larger number of iterations is used. Also, the error for the fine mesh converges to the same value as for the coarse mesh, provided the same design perturbation is applied. However, for the fine mesh it requires more iterations than for the coarse mesh to obtain convergence to the same error value. By means of an appropriate scaling of the axes of Figure 2.7 the curves in this graph can be made similar. For this purpose the quantity ∆x/(ηh2 ) is plotted along the horizontal axis. The error as defined above is divided by η and is plotted along the vertical axis. The result of this scaling is presented in Figure 2.8. From this graph it can be seen that at a value of 102 of the quantity ∆x/(ηh2 ) the error decreases sharply. For values of ∆x/(ηh2 ) smaller than 101 the error remains approximately constant. From this observation a stop criterion for the Laplacian smoothing can be proposed. This criterion could be formulated as ∆x ≤ ch2 η where c is a user-defined constant.. (2.26).

(27) Accuracy improvement of semi-analytical design sensitivities. 0. 19. η = 10−2 , coarse mesh η = 10−4 , coarse mesh η = 10−2 , fine mesh η = 10−4 , fine mesh. -0.5 -1. log ε. -1.5 -2 -2.5 -3 -3.5 -4 -4.5 -6. -5. -4. -3. -2. -1. 0. log (∆x/η). Figure 2.7: Error in first order refined semi-analytical design sensitivities using a coarse and a fine mesh and an initial relative perturbation of 10−2 and 10−4 .. 2.5. η = 10−2 , coarse mesh η = 10−4 , coarse mesh η = 10−2 , fine mesh η = 10−4 , fine mesh. log (ε/η). 2. 1.5. 1. 0.5. 0. -0.5 -4. -2. 0. 2. log (∆x/(ηh2 )). 4. 6. Figure 2.8: Error in first order refined semi-analytical design sensitivities using a coarse and a fine mesh and an initial relative perturbation of 10−2 and 10−4 ..

(28) 20. Chapter 2 0. log (d h2 /∆x). -0.5 -1 -1.5 -2 -2.5 -3 -3.5 -18. 10 elements 100 elements 1000 elements -16. -14. -12. -10. -8. log ∆x. -6. -4. -2. 0. Figure 2.9: Convergence of the distortions and perturbations during the Laplacian smoothing of a mesh with 10, 100 and 1000 elements.. It is somewhat remarkable that the horizontal axis has to be scaled using ∆x/(ηh2 ) and not using ∆x/η. The fact that 1/h2 is involved in this scaling causes that for accurate results many iterations are required for a fine mesh as ∆x must be iterated to a lower value as compared to a coarse mesh. In order to investigate this aspect in more detail, a one dimensional problem is investigated. For this problem the Laplacian smoothing is used. It turns out, that element distortions caused by the Laplacian smoothing converge significantly slower than the nodal perturbations. Moreover it follows that for larger numbers of iterations the convergence of the element distortions introduced by the Laplacian smoothing is proportional to ∆x/h2 which is illustrated in Figure 2.9. In this figure d h2 /∆x is plotted along the vertical axis and ∆x is plotted along the horizontal axis, where d represents the element “distortion”. Here the distortion denotes the largest relative change of the element lengths. It is seen, that the ratio d h2 /∆x becomes roughly constant for small ∆x. This observation indicates that the convergence of the design sensitivities is not governed by ∆x, i.e. the convergence of the nodal perturbations, but is determined by the convergence of the element distortions. As the element distortions turn out to converge much slower than the nodal perturbations (ratio is h−2 ), this causes a major handicap for the Laplacian smoothing. Results for the semi-analytical method are presented in Figure 2.10. Again the error decreases if the number of iterations is increased. However, results for the fine and the coarse mesh do not converge to the same error value. Also for the fine mesh it requires more iterations than for the coarse mesh to obtain convergence. By means of scaling of the axes of Figure 2.10 the curves in this graph can be made similar. Along the horizontal axis the same scaling as before has been used. The error is scaled by h2 /η and is plotted along the vertical axis. The result of.

(29) Accuracy improvement of semi-analytical design sensitivities 6. 21. η = 10−2 , coarse mesh η = 10−4 , coarse mesh η = 10−2 , fine mesh η = 10−4 , fine mesh. 5 4. log ε. 3 2 1 0 -1 -2 -6. -5. -4. -3. -2. -1. 0. log (∆x/η). Figure 2.10: Error in first order semi-analytical design sensitivities using a coarse and a fine mesh and an initial relative perturbation of 10−2 and 10−4 .. this scaling is presented in Figure 2.11. The use of h2 in the scaling of the error is in agreement with the analytical beam example from which it followed that the error is dependent on the number of elements squared.. 2.7.2 Cylindrical panel As a second example a cylindrical panel has been studied. The radius of the panel is 100 mm, the span (φ) is π/2 and the thickness is 0.6 mm. Young’s modulus is 3 × 105 N/mm2 and Poisson’s ratio is 0.3. The strip is clamped at one side and loaded by a distributed bending moment at the other side. The span of the panel (φ) is chosen as design variable. Design sensitivities have been computed for the tip displacement component perpendicular to the shell surface. The finite element mesh of the panel, the boundary conditions and the loading are shown in Figure 2.12. Along the curved sides symmetry boundary conditions have been applied. The analytical solution for the tip displacement is given by MR2 1 − ν2 u = 12 (φ cos φ − sin φ) . (2.27) Eh3 The first order design sensitivity follows as MR2 1 − ν2 du = −12 (φ sin φ) . (2.28) dφ Eh3 Again, logarithmic design sensitivities will be used, defined by u0l =. φ du = π2 /4. u dφ. (2.29).

(30) 22. Chapter 2. 3.5. η = 10−2 , coarse mesh η = 10−4 , coarse mesh η = 10−2 , fine mesh η = 10−4 , fine mesh. 3. log (εh2 /η). 2.5 2 1.5 1 0.5 0 -0.5 -4. -2. 0. 2. log (∆x/(ηh2 )). 4. 6. Figure 2.11: Error in first order semi-analytical design sensitivities using a coarse and a fine mesh and an initial relative perturbation of 10−2 and 10−4 .. Figure 2.12: Mesh of the cylindrical panel, boundary conditions and loading..

(31) Accuracy improvement of semi-analytical design sensitivities. 23. logarithmic 1st order design sensitivity. 3. 2.8. 2.6. 2.4. 2.2. 2 -11. SA RSA SA-Laplace RSA-Laplace -10. -9. -8. -7. -6. -5. -4. -3. -2. log relative perturbation. Figure 2.13: First order design sensitivities computed by means of the semianalytical method and the refined semi-analytical method, with and without Laplacian smoothing.. Again, first order semi-analytical and refined semi-analytical design sensitivities have been computed corresponding to different values of the relative perturbation (Figure 2.13). As in the previous example, it can be seen that if Laplacian smoothing is used, the range of relative perturbations resulting in accurate sensitivities is extended for both the semi-analytical and the refined semi-analytical method. A carefull inspection of the results depicted in Figure 2.13 shows that the results with and without smoothing converge to a different solution. The reason for this is very likely the fact that the response function is very sensitive with respect to the mesh layout and the mesh density. This is a consequence of the nearly inextensional bending mode and the small wall thickness studied here. To check this argument, response curves have been studied in the neighbourhood of the design point using both a boundary-layer perturbation and a smoothed perturbation. This study indeed confirms that the corresponding response curves differ, but intersect for corresponding meshes at the design point (Figure 2.14). In this figure the response functions are normalized by means of the analytical solution (2.27). For this problem the effect of the number of iterations on the accuracy of the design sensitivities has been studied. Both a coarse mesh (h = 0.064) and a fine mesh (h = 0.013) have been used. The error in the sensitivities has been computed for initial relative perturbations of 10−2 and 10−4 . This error has been defined as a relative error in which the numerical sensitivity is compared with the best numerical solution for the coarse and the fine mesh. Results for the semianalytical and refined semi-analytical sensitivities are presented in Figure 2.15 and Figure 2.16, respectively. The logarithm of the error has been plotted along the.

(32) 24. Chapter 2 1.25 1.2. response function. 1.15 1.1 1.05 1 0.95 0.9 Boundary (coarse mesh) Laplace (coarse mesh) Boundary (fine mesh) Laplace (fine mesh) Exact solution. 0.85 0.8 0.75 0. 0.02. 0.04. 0.06. 0.08. 0.1. relative perturbation. Figure 2.14: Response curves for the cylindrical panel using a boundary node approach and Laplacian smoothing for both a coarse and a fine mesh.. vertical axis. The logarithm of the maximum relative perturbation that occurred during Laplacian smoothing and scaled by the relative mesh size squared (h2 ) is plotted along the horizontal axis. For semi-analytical design sensitivities the error for the coarse mesh converges to the same value as for the fine mesh, provided the same design perturbation is applied. Refined semi-analytical design sensitivities show a clear convergence in case of a small design perturbation. Application of a similar scaling as for the strip example does not lead to a set of curves which become nearly identical. This must probably also be attributed to the fact that there is such a large influence of the model used on the accuracy.. 2.8. Discussion. Semi-analytical design sensitivities of finite element responses often suffer from severe accuracy problems in case of shape design variables. In the present paper it has been shown that the method which is used to create the perturbed mesh significantly influences the accuracy. This has also been found by Lindby and Santos (1997) who studied the effect of a boundary layer, boundary displacement and a mesh velocity method on the accuracy of the design sensitivities. By means of an analytical beam example it could be demonstrated that perturbation of all elements of the mesh results in more accurate sensitivities compared to a single perturbation of the tip element. This observation has been used to extend the approach to three-dimensional shell problems. By means of Laplacian smoothing the design perturbation can be distributed more uniformly over the domain which has a positive effect on the accuracy of the semi-analytical design sensitivities. The accuracy improvement of the semi-analytical design sensitivities has been.

(33) Accuracy improvement of semi-analytical design sensitivities. 25. 4 3.5. log ε. 3. η = 10−2 , coarse mesh η = 10−4 , coarse mesh η = 10−2 , fine mesh η = 10−4 , fine mesh. 2.5 2 1.5 1 0.5 -3. -2. -1. 0. 1. log (∆x/(ηh2 )). 2. 3. 4. Figure 2.15: Error in first order semi-analytical design sensitivities using a coarse and a fine mesh and an initial relative perturbation of 10−2 and 10−4 .. 0 -1. log ε. -2 -3 -4 -5. η = 10−2 , coarse mesh η = 10−4 , coarse mesh η = 10−2 , fine mesh η = 10−4 , fine mesh. -6 -7 -3. -2. -1. 0. 1. log (∆x/(ηh2 )). 2. 3. 4. Figure 2.16: Error in first order refined semi-analytical design sensitivities using a coarse and a fine mesh and an initial relative perturbation of 10−2 and 10−4 ..

(34) 26. Chapter 2. shown by means of several numerical examples. By using Laplacian smoothing the range of relative perturbations resulting in accurate sensitivities is extended for both the semi-analytical and refined semi-analytical method. If Laplacian smoothing in combination with a remainder is used the accuracy is improved considerably. A disadvantage is that the number of iteration steps may be large if full advantage is to be taken of the Laplacian smoothing. This particularly holds true if a fine mesh is used. The reason is that the convergence seems not directly related to the local mesh perturbations but to the element distortions. Implementation in an existing finite element program is simple. As the Laplacian smoothing is used at the end of the mesh creation procedure it can be simply added to an existing preprocessor. The implementation of the algorithm itself is also simple. The coordinates of the nodes of the mesh, the connectivity information and the number of boundary nodes should be available. A loop over all internal nodes of the mesh is performed and for each node the coordinates of the adjacent nodes are averaged and attributed to the current node. In the present paper it has been shown that the accuracy of semi-analytical design sensitivities can be improved by using Laplacian smoothing instead of a boundary node approach. However, this method has not been based on minimization of the error in the sensitivities. It might be possible that other mesh perturbation methods reduce the error even further. This will be the topic of future research..

(35) Reproduced from: K. Vervenne and F. van Keulen (2002), An alternative approach to response surface building using gradient information, 43rd Structures, Structural Dynamics, and Materials Conference, 22-25 April, Denver, Colorado, USA, AIAA-2002-1584. Chapter 3 An Alternative Approach to Response Surface Building using Gradient Information Abstract Response surface techniques are widely used to fit experimental data. If the data is obtained from numerical analyses, then gradient information may be available inexpensively. This sensitivity data could be used to improve the accuracy of the response surface or to reduce the number of points required. The gradient information can be incorporated using (iteratively) weighted least squares. The present paper presents an alternative approach consisting of two steps. In the first step the function and gradient data is fitted separately, resulting in a response surface for both the function itself and one for every derivative. In the second step these response surfaces are combined to form a single response surface. It is shown that including gradient information in the response surface may be considered as a multi-objective optimization problem, in which a trade-off between accuracy in function values and gradients has to be made explicitly. The alternative approach will be demonstrated and studied by means of a series of simple examples. Key words: response surface, fitting, regression, design sensitivities, gradients, least squares.. 3.1. Introduction. Response surfaces are widely used to fit data obtained from experiments. In literature a large number of well-crystallized methods can be found to build response 27.

(36) 28. Chapter 3. surfaces for this purpose (Myers and Montgomery 1995; Khuri and Cornell 1996). During the response surface building, an approximation function is fitted to given experimental data. The parameters of the approximation functions are estimated by means of a least squares approach. In the present paper the building of response surfaces will be discussed based on data obtained from numerical experiments. These response surfaces might be used for optimization procedures, function estimation or robust design. The response surface approach is particularly attractive when the function evaluations are expensive in terms of computational cost. Another advantage of response surface approximations is that possible noise in function evaluations is smoothened. An important requirement of a response surface is its capability to represent the data accurately. The accuracy of the response surface might be improved by using derivative information, i.e. derivatives of the response functions with respect to the (design) variables. Especially in cases where derivative information is available at low computational cost, use of this information might be advantageous. For example in case of finite element analyses, derivative information is often available inexpensively. By means of the semi-analytical method the sensitivities can be obtained at a fraction of the time of the analysis itself (Haftka 1993). In literature several methods to include derivative information in response surface approximations can be found. In (van Keulen, Liu, and Haftka 2000; Lauridsen, Vitali, van Keulen, Haftka, and Madsen 2001) sensitivities are included into the response surface by an enhanced weighted least squares formulation. A problem associated with this method is the selection of weights used to indicate the relative importance of the sensitivities compared to the function values. In (van Keulen, Liu, and Haftka 2000) several weighting schemes are discussed. The simplest weighting method is standard normalization of both function values and derivatives. An alternative scaling method uses logarithmic design sensitivities. Finally, the weights might be selected by an iterative algorithm in which the weights are based on an estimation of the covariance matrix (Lauridsen, Vitali, van Keulen, Haftka, and Madsen 2001). In (Malkov and Toropov 1991) design sensitivity information is also used to obtain an approximate model of system behaviour. The fitting is based on a weighted least squares method. The approximate models are used in an iterative procedure for solving structural optimization problems. Here the weights are based on normalization. Particularly the norm of the gradient has been used. In (Rijpkema, Etman, and Schoofs 2000) sensitivities are used in response surface building by a weighted least squares approach. The weights are determined iteratively by estimation of the covariance matrix. Furthermore, kriging models are extended by derivative information. The authors conclude that the efficiency of the model training can be improved by using sensitivity information. However, if the response is discontinuous they recommend not to use derivative information..

(37) Response surface building using gradient information. 29. If the derivative information is obtained by global finite difference techniques it is also not advantageous to use sensitivities. Here it is assumed that calculation of the sensitivities requires a full re-analysis. Instead, response data should be obtained for a larger set of design points. The method studied here includes derivatives into the response surface by means of a two-step approach. In the first step the function data and derivative data are fitted by separate response surfaces. In the second step a single response surface is constructed based on the response surfaces of function and derivatives. This response surface might be of higher order than the individual response surfaces for the function and gradients. It will be demonstrated that the combination of function and derivative data in a single response surface can be considered as a multi-objective optimization problem. If the response surface building is treated like this, a good compromise between accuracy in function values and derivatives can be made explicitly. The outline of the present paper is as follows. An alternative approach to response surface building and the multi-objective formulation are described in Section 2. Examples are studied in Section 3. Results are discussed in Section 4.. 3.2. Response surface building. In the present section derivatives enhanced weighted least squares will be summarized. After that an alternative (two-step) approach for response surface building using sensitivities will be discussed. Finally, the mixing of function values and derivatives will be formulated as a multi-objective optimization problem.. Derivatives enhanced weighted least squares In this section derivatives enhanced weighted least squares will be summarized briefly. A single response function y(x) will be considered, where x ∈ Rn . It will be assumed that function value and derivative data is available in s points xi (i = 1 . . . s). The function values will be denoted by yi0 = y(xi ) i = 1 . . . s,. (3.1). and the derivatives will be denoted by yi j =. ∂y (xi ) ∂x j. i = 1 . . . s, j = 1 . . . n.. (3.2). A compact notation is introduced by yTi = [yi0 . . . yin ] .. (3.3).

(38) 30. Chapter 3. The data in the s plan points will be approximated by a response surface which can be written as yˆ i = Ai β. (3.4) Here the vector β represents the coefficients of the response surface. The matrix Ai contains the approximation functions and its derivatives with respect to the coordinates xk (see (van Keulen, Liu, and Haftka 2000) for more details). The coefficients will be determined using a weighted least squared approach. In contrast to standard weighted least squares, information on both function values and derivatives will be used. In order to carry out the least squares, an error function is defined by s. E = ∑ (yi − yˆ i )T Wi (yi − yˆ i ).. (3.5). i=1. Here, Wi is a weighting matrix which can be used to attribute different weights to function values and derivatives. Within the present paper it is assumed that Wi is identical for all points. If no weighting is used Wi is equal to I. The coefficients are computed by minimizing the error function. This results in a linear set of equations which can be solved for the coefficients β.. Two-step approach This section describes an alternative approach to response surface building using function evaluations and derivative information. It will be assumed that for a given plan of experiments function values and design sensitivities are available. The method described here can be viewed as a two-step approach. In the first step the following response surfaces are constructed. One based on the function values only and one for each derivative. Hence n + 1 intermediate response surfaces are constructed. These response surfaces can be build using existing methods for response surface building. In the second step the response surfaces for the function values and the derivatives are combined to form a single response surface. This is done in such a way that the final response surface fits as good as possible to both the function value response surface and the derivative response surfaces. One way to achieve this is to minimize the difference between the final response surface and the response surface of function values and at the same time minimizing the differences between the derivatives of the final response surface and the response surfaces for the derivatives. This minimization procedure is performed by using a least squares approach. In the sequel the response surface of the function values will be denoted by fˆ0 and the response surfaces of the derivatives by fî . The final response surface will be denoted by pˆ0 and its derivatives by pî . In order to carry out the minimization.

(39) Response surface building using gradient information. 31. procedure described above, an objective function is formulated which is defined by Z n 2 E = ∑ wi pî − fî dV, (3.6) i=0. V. where wi are weight factors. The integral is evaluated over the total design space which is spanned by the n (design) variables. It is noted that the objective function given above is the most simple one; the weight factors might be a function of the design variables and thus be placed inside the integral. It will be assumed that the final response surface and its derivatives depend linearly on the coefficients: pî (x, b) = pˆ Ti (x)b.. (3.7). Here x represents the design variables and b is a set of coefficients. By means of a least squares fitting procedure, the coefficients of the final response surface are determined. Substitution of (3.7) in (3.6) and differentiating with respect to the coefficients results in T n ∂E = ∑ 2wi [Pi b − ai ] , (3.8) ∂b i=0 Z. with Pi =. V. pˆ i pˆ Ti dV. (3.9). pˆ i fî dV.. (3.10). Z. and ai =. V. To minimize the objective function, the derivative of the objective function with respect to the coefficients is equated to zero which gives n. ∑ 2wi [Pib − ai] = 0.. (3.11). i=0. Equation (3.11) is a linear set of equations which is solved for the coefficients of the final response surface. An advantage of the method described here is that in the first step of the procedure the function values and derivatives are fitted separately. Advanced techniques for building response surfaces based on function values only have been developed. These techniques can be used to build the response surfaces in the first step. In the second step of the procedure the two types of response surfaces are mixed. The selection of weights that determines the relative importance of the derivatives with respect to the function values is not obvious. A possible solution to this problem is to consider the second step as a multi-objective optimization problem. This will be the topic of the next section..

(40) 32. Chapter 3. E1. E0. Figure 3.1: Pareto optimal set for a single design variable, E0 is the error in the function value and E1 is the error in the gradient.. Multi-objective formulation Construction of a response surface using both function values and derivatives can be considered a multi-objective optimization problem. The objective functions which are to be minimized simultaneously are the errors in the function value and in the derivatives. This results in n + 1 objective functions, where n denotes the number of design variables. Often, these objectives are competing; improvement of one objective leads to deterioration of another. Using the notation introduced previously, the objective functions are given by Z 2 Ei = pî − fî dV, (3.12) V. where E0 denotes the error in the function value and Ei , i = 1 . . . n denote the errors in the gradients. In equation (3.6) the objective functions are combined into a single objective function using weights corresponding to each objective function. For a particular choice of the weights a trade-off between errors in function values and derivatives is obtained. By variation of the weights used in the regression, a Pareto optimal set can be generated. This has been illustrated in Figure 3.1 for the case of one design variable. In order to make a justified choice for the weights several criteria might be used. For example, the sum of squares of the errors could be minimized with respect to the weights. In that case the objective function to be minimized is given by n. ENorm = ∑ Ei2 . i=0. (3.13).

(41) Response surface building using gradient information x -1 1. f (x) 0 4. 33. d f /dx(x) 2 6. Table 3.1: Function and derivative data in the two plan points x = −1 and x = 1.. Note that the errors Ei can only be computed after solving (3.11) for the coefficients and substitution in (3.12). Alternatively, the maximum of the errors Ei could be minimized, which results in the following objective function EMax = max. Ei .. (3.14). Scaling of the errors Ei might be required because their order of magnitude may differ. In (van Keulen and Vervenne 2002) a scaling method is proposed. The scaling of the function values is defined by min 2 E0 /(ymax 0 − y0 ) ,. (3.15). and the derivatives are scaled by min 2 4Ei /(|ymax i | + |yi |) ,. i 6= 0.. (3.16). and ymax Here ymin i i ( j = 0 . . . n) are the minimum and maximum function values and gradients.. 3.3. Examples. The present section presents two examples of response surface building using gradient information. By means of these examples the alternative approach is compared with a direct approach. Also the effect of the applied multi-objective criterion is studied. Finally, for the first example the influence of the number of points is examined.. 3.3.1 Example I In this section the method described in the previous section will be applied to the fitting of a polynomial in one dimension. The function to be approximated is a third order polynomial with coefficients equal to 1, defined on the domain [−1; 1]: f (x) = 1 + x + x2 + x3 .. (3.17).

(42) 34. Chapter 3 4 3.5. c1. coefficients. 3 2.5. c0. 2 1.5 1. c2. 0.5 0 0. 0.2. 0.4. 0.6. 0.8. 1. w. Figure 3.2: Coefficients c0 , c1 and c2 of the final response surface as function of the weight factor w which determines the relative importance of function value compared to derivative.. It will be assumed that the function values and derivatives are known for the plan points x = −1 and x = 1. These are listed in Table 3.1. As a first step the function and derivative data is approximated by linear functions. For the function values this results in fˆ0 (x) = 2 + 2x. (3.18). and for the derivatives the linear approximation is given by fˆ1 (x) = 4 + 2x.. (3.19). As a second step these linear approximations are combined into a single response surface. A second order polynomial will be used for this response surface: p(x) ˆ = c0 + c1 x + c2 x2 .. (3.20). By means of a least squares approach the coefficients of this polynomial are determined. The error function that is minimized is E = wE0 + (1 − w)E1 ,. (3.21). where E is the total error and w is a weight factor. The error in the function value, E0 is defined by Z 1 2 fˆ0 − pˆ dx, E0 = (3.22) −1. and the error in the derivative, E1 is defined by Z 1. E1 =. −1. 2 fˆ1 − d p/dx ˆ dx.. (3.23).