Robust optimization strategies for sheet metal springback compensation

Robust optimization strategies for sheet metal springback compensation XIII International Conference on Computational Plasticity. Fundamentals and Applications COMPLAS XIII E. O˜nate, D.R.J. Owen, D. Peric & M. Chiumenti (Eds) ROBUST OPTIMIZATION STRATEGIES FOR SHEET METAL SPRINGBACK COMPENSATION A. Maia∗, E. Ferreira†, M.C. Oliveira†, L.F. Menezes†and A. Andrade-Campos∗ ∗Department of Mechanical Engineering, Centre for Mechanical Technology & Automation, GRIDS Research Group, University of Aveiro Campus Universit´ario de Santiago, 3810-193 Aveiro, Portugal e-mail: [email protected], web page: http://www.ua.pt/ †Center for Mechanical Engineering of the University of Coimbra (CEMUC) Department of Mechanical Engineering, University of Coimbra P´olo II, Rua Lu´ıs Reis Santos, Pinhal de Marrocos, 3030-788 Coimbra, Portugal e-mail: cem[email protected], web page: http://www.uc.pt/en/iii/research centers/CEMUC Key words: Sheet metal forming, Springback Compensation, Finite Element Model Updating Strategy, Robust Optimization strategies, NURBS parametrization Abstract. Sheet metal forming is a major industrial process, mainly due to its cost efficiency after the establishment of the process design. However, the process design from tools geometry to load conditions is not straightforward, as a consequence of the side effects associated with sheet metal forming. The emphasis in this area goes to the springback effect or elastic recovery, which is one of the main causes of part’s inaccuracy, demanding tool compensation. This work proposes to compare different robust optimization strategies to sheet metal forming springback compensation. The methodology adopted resorts to Response Surface Method (RSM), as well as to Finite Element Model Updating (FEMU) strategies, to adjust the design variables. These include the tools’ surfaces, which are parametrised with NURBS. These strategies are then compared using the U-Rail benchmark. The results achieved reveal a reduction of 99% on the geometrical error of the final piece for the best methodology. 1 INTRODUCTION The evolution observed in the industry, where new pieces are constantly needed, creates the necessity to accelerate the design stage of new forming processes while maintaining its accuracy. Concerning springback compensation, this means that the necessary tools adjustments and other variables need to be swift and rigorous. This challenge is not efficiently answered by the traditional ”trial-and-error” process. Therefore, the employment of optimization strategies, among which statistical methods, such as the Response Surface Method (RSM), or Finite Element Model Updating (FEMU) strategies associated 783 A. Maia, E. Ferreira, M.C. Oliveira, L.F. Menezes and A. Andrade-Campos to optimization algorithms, have been studied. Though there are commercial software [1] that can predict the final geometry of a formed sheet given the initial variables and adjust them to achieve the desired piece (direct and inverse problem respectively), their efficiency is still not optimal. The work hereby presented proposes to implement an integrated methodology that searches for the set of design variables that better compensates the springback effect. The results obtained are then checked for feasibility and a comparison between the proposed optimization strategies is made. The design variables considered are the control points of the NURBS that define the forming tools geometry and the Blank Holder Force (BHF). 2 SPRINGBACK COMPENSATION STRATEGIES A single-step sheet metal forming process is divided in three phases: (i) the initial positioning of the tools; (ii) the forming by the tools and (iii) the removal of the tools with subsequent springback. Springback or elastic recover is an undesired side-effect of sheet metal forming which is proportional to the ratio between residual stresses and Young Modulus [2]. This relation leads to a first approach for the minimization of this side-effect: an increase on the plastic deformation, done by incrementing the BHF and/or other restraining forces [3]. However, this approach may induce serious quality problems in the final part as this leads to a reduction on the material flow during the forming phase. This may generate necking of the blank and consequent fracture. Another possible approach are Springback Compensation methodologies. These methodologies acknowledge the springback as unavoidable. Thus instead of trying to reduce it, they use the springback so that the final piece has the desired shape. This type of approach comprises several strategies. The Displacement Adjustment method (DA) [4] consists on the displacement of the tool’s surface on the opposite direction of the springback, which can be done in just one step [3, 5, 6] or iteratively [7]. The later concept has revealed itself more effective in areas strongly affected by plastic deformation. In spite of its good practical results, this methodology has some drawbacks, such as the difficulty in aligning the CAD model and the forming piece or the tendency of the compensated surface to become rougher [3]. The former problem may be solved by the Smooth Displacement Adjustment (SDA). This methodology approximates the geometric error by smooth functions, using boundary constraints if restrictions are needed [5]. Other improvements on the DA strategy is the Comprehensive Compensation (CC), where instead of adopting the opposite direction of the springback, the best direction and magnitude of the adjustment is computed. This approach combined with the iterative DA has achieved promising results [6]. The Direct Curvature Method (DCA) is another possible approach. Instead of treating the tool as an whole, it changes different parts of the tool independently, through a dynamic compensation factor [3]. The computation of the correct compensation factors / adjustments can be done through a ”trial-and-error” methodology using numerical trials to achieve the optimal adjustment, or using optimization algorithms. 784 A. Maia, E. Ferreira, M.C. Oliveira, L.F. Menezes and A. Andrade-Campos Optimization algorithms are methods to compute the set of input variable that minimize or maximize a cost-function. These algorithms may be classified in (i) Natureinspired algorithms, (ii) Gradient-based methods and (iii) Artificial intelligence (AI). While the Gradient-based algorithms present a faster convergence and are fairly efficient, they may be trapped in local minima. The Nature-inspired methods statistically converge for the global minimum requiring however many evaluations of the cost function, which leads to high computational costs. The AI algorithms also have the disadvantage of requiring massive quantities of data for their training. However, when well calibrated, they provide very accurate results. In sum, all methods have advantages and disadvantages, being the choice of their employment strongly dependent on the specific application desired and on the amount of data available. There is also the possibility to combine several strategies and approaches in order to take advantage of their best features [8]. In this work, a similar approach to the DCA is adopted. However, the tools modelling and adjustments is done through NURBS, being the design variables the NURBS control points. The geometric evaluation however remains using geometrical parameters, as these better express the desired shape for the component. 3 METHODOLOGY AND IMPLEMENTATION In order to find the design variables that optimize (minimize) the cost function, the process presented in Fig. 1 is followed. The cycle starts by pre-processing the original input files which store the initial design variables, including the desired part design and the initial tool design, defined by IGES and Mesh files (one per tool). Following, the simulation software outputs the numerical results corresponding to the formed metal sheet. An evaluation of these results is then computed, based on the geometrical error between the obtained piece and the reference one. This cost function is subsequently introduced into the optimization algorithm, which computes the new design variables. At this stage, an evaluation of the whole process takes place: if the convergence criteria of the algorithm is met, then these variables are outputted as optimised values; otherwise the input files are updated with the new variables proposed, pre-processed and input on the simulation software, starting a new cycle. Figure 1: Process Schematics 785 A. Maia, E. Ferreira, M.C. Oliveira, L.F. Menezes and A. Andrade-Campos It should be noticed that this cycle only applies to the Direct Search and Gradientbased algorithms. For the RSM, first all the FEM-simulations are run for a predetermined set of variables combinations, and only at the end the optimization is applied. Each action listed is performed by an individual block of code. This option for modularity lends flexibility to the code and eases the error detection and the code evolution. 3.1 Finite Element Analysis (FEA) The simulations are carried out using the finite element code DD3IMP [9], developed specifically to simulate sheet metal forming processes, using an updated Lagrangian formulation and a predict-correction scheme to determinate the equilibrium state [10]. A problem inherent to the use of the FEM is the numerical noise due to round off errors, mesh discretisation and instability of contact conditions [11]. This noise may lead to different meshes, variations on the number of increments of the simulations and, consequently, to different results and evolutions of the optimization algorithm. In extreme cases, the noise can create local minima trapping the optimisation process, and lead to inaccurate optimization results. As the noise creates more local minima of the objective function, it is necessary to account for its impact in the optimization evolution. To this end, a multi-start strategy in conjunction with the chosen optimization algorithms is employed. 3.2 Springback compensation through tool design This implementation relies on Finite Element Model Updating (FEMU) strategy [12]. Once started, the optimization algorithm does its iterations resorting to a simulation every time a new value for a set of design variables is necessary. To do so, the algorithm outputs a set of design variables and automatically initiates a cycle of rewriting and pre-processing new input files, performing a simulation and evaluating the obtained piece, using this value to the next iteration. Although this approach is computationally expensive and subject to numerical noise, it is reliable in the sense that it is not based on interpolations of the cost function behaviour. Instead, it computes the actual values of the measure variables for each set of design variables. In the scope of this article a least-squares Gradient-based algorithm, a direct search algorithm and a linear Response Surface Method are used. The results achieved by all the methods are computed and discussed. It should be noticed that all the design variables are considered as relevant inputs of both the interpolation and the optimization algorithms. 3.2.1 Linear Response Surface Method (RSM) This approach resorts to an interpolation of a set of pre-existent residues in order to estimate equations that describes the behaviour of the cost function. Previously to the optimization process, a Sensitivity Analysis was conducted to ascertain the general behaviour of the blank upon perturbations on the process design variables and to help to choose the optimization strategies starting points. Having all the data, a simple Linear 786 A. Maia, E. Ferreira, M.C. Oliveira, L.F. Menezes and A. Andrade-Campos Fitting is applied, resulting on a set of linear equations that describe a planar Response Surface of the geometric evaluation parameters. 3.2.2 Direct Search Algorithm The Nelder-Mead Simplex Algorithm [13], a direct search method, relies on the construction of a simplex of N+1 vertices, for a N-dimensional problem. Then, it iteratively replaces its vertices for new ones with lower values of the cost function. Its main advantage is independence of the gradient of the cost function or any approximation. To study the influence of the unconstrained nature of the Nelder-Mead algorithm a simulation where there is no constraints to the values of the input variables, and another where a variable transformation is applied are performed. The transformation is applied between the unconstrained variable xi, to be optimized, and the constrained variable Xi, subjected to a maximum Xmax iand a minimum Xmin i, such as Xi=Xmin i+(Xmax i−Xmin i)exp(xi) exp(xi)+ exp(−xi).(1) This transformation is applied before each simulation, so that it is performed with feasible variables. It is also necessary to apply the inverse transformation at the beginning of the optimization, meaning that the constrained initial design variables Xioshould be transformed in the unconstrained ones xio. The tolerance of 10−5is employed, along with the standard values to the algorithm coefficients: α=1., γ=2. and β=δ=0.5. 3.2.3 Least-Squares Gradient Based Algorithm The Levenberg-Marquardt is a gradient-based method, similar to the Newton-Raphson, having however a stabilization parameter, µk, in order to improve the algorithm’s behaviour around minima. As this particular problem doesn’t have an a-priori gradient formula, the Jacobian matrix is computed at each iteration through forward finite-difference calculations. A major difficulty on this methodology for the sensitivity compensation is its sensitivity to noise. This noise not only affects the evaluations of the cost functions (possibly even introducing extra local minima), but also the Jacobian construction, which is crucial for gradient-based methods. In order to overcome this obstacle, a multi start strategy is implemented, equivalent to the one described on the Nelder-Mead subsection. The influence of the noise on the Jacobian, however, needed another approach. An increase on the finite-differences perturbation is made so that, even affected by noise, the difference translates the real trend of the cost function, which overcomes the noise. The chosen step for this particular application is 2% of the variable value. This algorithm implementation has a tolerance of 10−5, identical to the one adopted by the Nelder-Mead. 787 A. Maia, E. Ferreira, M.C. Oliveira, L.F. Menezes and A. Andrade-Campos 3.3 Process Parametrization Classical types of representation, such as B´ezier curves or geometrical parametrizations (lengths, radius and angles), though being more intuitive to the human user, do not have the necessary flexibility to parametrize very complex geometries or became very heavy to compute. More advanced parametrizations such as the T-splines, though presenting important advantages, are still not spread in industry and present few in-use software. The use of NURBS [14] as parametrization presents a good balance. Their flexibility and accuracy on representing complex shapes have lead to their adoption and standardization at an industrial level, being used on a wide range of applications and software. In this work, NURBS are used as parametrization methodologies to define the design variables: their control polygon coordinates and weights. All the data concerning NURBS control points is fully codified in standard IGES files, widely applied in industry. The IGES files pre-processing occurs in two phases (i) different IGES files with the new design variables are written; (ii) the constraints imposed by the problem under analysis are verified and the necessary meshes generated, using GiD software through a batch file. The Blank Holder Force (BHF) is an additional design variable which has to be closely monitored due to possible necking. This is relevant in order to ascertain the feasibility of the results found at the blank structural level. 3.4 Cost Function Formulation The main goal is to minimize the geometrical deviation between a formed piece and the reference one. This deviation may be measured by means of the euclidean distance between corresponding pairs of nodes of the formed blank. Alternatively it can measure the distance between those nodes along one well defined direction (e.g. X-axis) [6] or rather as the error between agreed sets of parameters (such as the difference between a certain angle). The use of distances requires a pairing algorithm between optimised piece and reference nodes, which can be costly. More often, the option of comparing a set of deduced geometric parameters is enough. In this work, the cost function is defined as E(x)=1 n n ∑ i=1 bi(xi−xiref)2,(2) where xiand xref iare the observed and desired values and biis the weight associated to the ith parameter difference. The weight is also used as a scale factor in order to level the differences of magnitude orders. 4 CASE STUDY: U-RAIL The springback compensation strategies described on the previous section are compared using the U-Rail benchmark. This particular case study is chosen due to large springback effects which, if not dully compensated, may cause serious quality problems. 788 A. Maia, E. Ferreira, M.C. Oliveira, L.F. Menezes and A. Andrade-Campos 4.1 Case Study Features Fig. 2 presents the U-RAIL problem, as described in [15], which is taken as reference to this work. The original blank has 300x300x0.8 [mm] but it can be simplified assuming plane strain conditions. Thus only half of the length and a strip of 10 mm width is simulated. This is done to reduce the computational time. The processes is implemented on a computer equipped with an Intel Core I7 Processor with 8 cores and with 8 GB RAM, where each simulation takes an average of 711 seconds to be computed. The blank material has the properties listed on Table 1. This table also shows the number of elements used in the FEM discretization. The tools are assumed as rigid and handled by means of Nagata patches [16]. A friction coefficient of 0.15 is used. (a) Tools and design variables (b) Perspective view of the reference piece (c) Desired piece geometry Figure 2: Schematic of the experiment relevant angles and points Property Value Young Modulus [MPa] 206629 Poisson’s Ratio 0.298 Swift Hardening Law (K [MPa], n, ϵ0)488.35, 0.24, 0.015 Hill48 (F,G,H,N) 0.63974, 0.60976, 0.39024, 1.43693 Number of elements 990 (length) x 1 (width) x 3 (thick.) Table 1: Material Properties and model parameters In order to verify if the conditions adopted on DD3IMP accurately capture the real behaviour of the metal sheet forming, the results are compared for the reference condition. The flange angle computed by DD3IMP is of 13.66o, which has an error of 0.07% to the values measured in [15]. The control points of the NURBS considered as design variables are showed in Fig. 2a: P1 to P4. Some constraints on their mobility are assured in order to guarantee some tools’ geometric features, such as orthogonality between vertical and horizontal surfaces or the fact that the concordance between those surfaces is described by a radius. They can be written as P3x= (30 −P1z); P4z= (28.85 −P2x),(3) 789 A. Maia, E. Ferreira, M.C. Oliveira, L.F. Menezes and A. Andrade-Campos where the subscript indicates the direction. Therefore only three design variables are studied: the coordinate zof point P1 to define the die shape, the coordinate xof P2 to define the punch geometry, and the Blank Holder Force. In cases where some limitation of the parameters is necessary to the optimization algorithm, the input space considered in [15] is adopted. The same is valid for the values adopted in Eq. 1 for Xmax iand Xmin i. 4.2 Case Study Evaluation (a) Points and angles (b) Zoom vertical surface Figure 3: Error measurements As the main objective is to ensure the orthogonality and planarity of the surfaces, instead of computing the euclidean distance between matching pairs of points on the two pieces, the geometrical error is computed resorting to the errors of the angles θ1and θ2 and the distance hto the respective reference values (see Fig. 3). In order to evaluate these variables, auxiliary points are defined in the component as shown in Fig. 3. Their coordinates are used to evaluate the deviations on the final piece trough α= tan (Ez−Dz Ex−Dx);θ1=180 −tan (Bz−Az Bx−Ax);h=√(Cx−Fx)2+ (Cz−Fz)2.(4) where the letter in subscript describe the coordinate of the point in case. The angle θ2is measured considering the trigonometric relationships between θ1and α,i.e., θ2=θ1−α. (5) The values of the angles θ1,θ2and distance h, which are evaluation variables, are used as inputs on Eq. 2 resulting in E(θ1,θ 2,h)=1 3((θ1−90)2+ (θ2−90)2+(h∗10)2).(6) As both angles are desired to be orthogonal, the reference values for both is 90o. The curvature of the vertical surface is supposed to be null, being the reference value of h 0 mm. Additionally, as the scale of the curvature values is a magnitude lower than the angle errors, this variable is multiplied by 10, so that its influence is not overlooked. 790 A. Maia, E. Ferreira, M.C. Oliveira, L.F. Menezes and A. Andrade-Campos 5 RESULTS DISCUSSION 5.1 Sensitivity and Noise Analysis Fig. 4 presents the sensitivity analysis performed to also evaluate the numerical noise. The least sensitive variable is the P2xthat defines the punch radius, what is in agreement with literature [15]. The other variables present a non negligible influence on the final geometry of the piece, being the BHF the one who influences the angles the most. Fig. 4c shows that an increase on the BHF values lead to an evolution of the geometric evaluation parameters of the piece to the desired references. The same is verified, though in a very smaller scale for low values of the P1z(small die radius). On this analysis, only his visibly affected by numerical noise. However, with these noise levels, it is expectable that the basic strategies employed in this work and presented in Subsection 3.2, are able to overcame it. In fact, the process presented in Fig. 1 is classified as robust, as it is able to overcome noise and failures on simulations. Those failures only occurred when the design variables fell into values that lead to impossible tools (e.g. tools with null negative radius). In those cases, the process follows to the next cycle with a new set of variables. Six starting points are used in this work. Due to its low sensitivity, the P2xstarting point value remained constant, being the points the combinations of the values −3, −5 and -7 for the P1zand 210 and 270 kN for the BHF. (a) P2x= 23.85mm; BHF = 210kN (b) P1z= -5.0mm; BHF = 210kN (c) P1z= -5.0mm; P2x= 23.85mm Figure 4: Sensitivity Analysis Trials: 5.2 Optimization Results The following expressions are obtained using the Linear RSM: θ1=−0.2014 ·P1z−0.0343 ·P2x−3.60 ×10−5·BHF + 102.65; (7a) θ2=−0.2076 ·P1z+0.0117 ·P2x+2.21 ×10−5·BHF + 81.00; (7b) h=−0.0129 ·P1z−0.0020 ·P2x−1.26 ×10−6·BHF +0.4502.(7c) The minimum of equation 6 using the previous definitions for θ1,θ2and hfall in a solution that is not feasible, considering that the punch and the die would intersect each 791