- Email: [email protected]
Structural design of aircraft skin stretch-forming die using topology optimization
Journal of Computational and Applied Mathematics 246 (2013) 278–288
Contents lists available at SciVerse ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
Structural design of aircraft skin stretch-forming die using topology optimization Ji-Hong Zhu a,∗ , Xiao-Jun Gu a , Wei-Hong Zhang a , Pierre Beckers b a
Engineering Simulation and Aerospace Computing (ESAC), Northwestern Polytechnical University, Xi’an, Shaanxi 710072, China
b
LTAS-Infographie, University of Liège, Liège 4000, Belgium
article
info
Article history: Received 23 January 2012 Received in revised form 29 August 2012 Keywords: Stretch-forming die Light-weight design Topology optimization Linear interpolation Casting constraints
abstract This paper demonstrates a topology optimization method for the large scale stretchforming die design. The implementation of the topology optimization depends on the surface loads obtained from the numerical simulation of the sheet metal stretch-forming process. Typically, the design variables are defined based on a linear interpolation to generate directional structures satisfying the manufacturing requirement. To further validate the design procedure, mathematical derivations have proved that the design sensitivities are strictly continuous. More explanations and numerical tests are presented to show the variation of the surface loads versus the stiffness of the die. A stretch-forming die design example is solved on account of gravity and surface loads. The final solution is compared with the traditional design. The numerical results have shown that the topology design can improve the stiffness and strength of the stretchforming die significantly. © 2012 Elsevier B.V. All rights reserved.
1. Introduction In the past two decades, the techniques of topology optimization were developed to generate reasonable structural configurations by re-distributing the material in the design domains with the prescribed loads and boundary conditions [1–4]. Till now, one of the most challenging works involved in the research of topology optimization is the industrial application, where complex practical load cases have to be taken into account and different manufacturing constraints have to be satisfied. The skin stretch forming die discussed in this paper is one of the key equipments of aircraft manufacturing. It basically consists of a certain stretch working surface of specific form and the supporting structure that interconnects the working surface for its integrity. Nowadays, the development of large cargo and passenger aircraft has brought a series of new challenges. On the one hand, the size and weight of the stretch-forming dies increase significantly. As the traditional die is always designed extremely stiff, it results in excessive costs as well as difficulties in transportation and installation. On the other hand, the structural design of the die must ensure the accuracy of the stretch-forming process. As a result, the light-weight design of the large scaled stretch-forming die using the structural optimization methods becomes necessarily important. In this paper, the topology optimization method is used to design the supporting structures on account of the gravity of the die itself and the surface loads obtained from the numerical stretch forming procedure.
∗
Corresponding author. E-mail address: [email protected] (J.-H. Zhu).
0377-0427/$ – see front matter © 2012 Elsevier B.V. All rights reserved. doi:10.1016/j.cam.2012.09.001
J.-H. Zhu et al. / Journal of Computational and Applied Mathematics 246 (2013) 278–288
279
(b) The pre-stretching.
(c) The die’s uprising.
(a) A stretch-forming system and its working procedure (www.custompartnet.com).
(d) The final stretching.
Fig. 1. An illustration of stretch-forming procedure.
The structural optimization design of the stretch-forming die proposed in this paper is actually carried out as a twostep way. First, the surface loads on the die have to be numerically predicted with the performing of non-linear finite element analysis. Some typical and extreme load cases will be selected for the topology optimization. Currently, most of the existing works on the stretch-forming were concentrated on the numerical simulation of different kinds of forming process to improve the forming precision [5]. For example, Yan and Klappka [6] studied the spring-back behaviors of panel forming productions using the multi-point stretch forming technique. Wang et al. [7] carried out a series of numerical simulations. The shape error and non-uniformity in thickness distribution of the sheet part were analyzed. Jiang et al. [8] combined the stretch forming with the roll forming procedure. The non-linear Finite Element simulation had shown more benefits that better forming precision and lower stress level were obtained. Till now, the development of the numerical simulation of the metal sheet stretch forming procedure has been satisfied. As a result, we choose to use the explicit time-integration method in the platform of ABAQUS to obtain the different surface load cases. Second, based on the contact forces, topology optimization will be implemented to design the supporting structure satisfying the stiffness of the die and the requirement of weight-saving. Currently, the application of topology optimization to the contact problems is not a new issue. For example, Petersson and Patriksson [9] proposed a subgradient algorithm to solve the topology optimization of sheets in contact. The convergence and the computational efficiency of the optimization were in-depth discussed. Recently, in the work of Strömberg and Klarbring [10], the unilaterally constraints were simulated with contact conditions. Both 2D and 3D problems were optimized and reasonable results were obtained. In the existing works, the structural optimizations of the metal forming dies are mostly found designing the stamping die. For example, Sheu and Yang [11] predicted the pressures on the working surface of a stamping die. The inner structures were then designed by modifying the size and layout of the stiffeners. The optimized results turned out to be completely different from the traditional design, i.e. a uniform distribution of the stiffeners. Nilsson and Birath [12] simulated the variation of the surface loads during the stamping process by means of time integration, where the process of lifting and stamping are taken into account. Topology optimization was then applied to save the weight by maintaining the structure rigidity and strength. Similarly, Xu and Tang [13] developed the topology optimization procedure for the inner structure design of the stamping die based on the platform of LS-Dyna and Hyperworks, where 28% of the total weight is saved. Based on the previous works, a stretch-forming die of a large cargo aircraft is designed in this paper. The stretch-forming process is first simulated based on ABAQUS to obtain the surface loads which will be defined as numbers of weighted load cases. By deriving the design sensitivities and applying the constraints on the material distribution and the total material cost, topology optimization of the inner structure is then carried out to find the optimized structural layout. 2. Basic formulation 2.1. Stretch-forming simulation A typical stretch-forming procedure as illustrated in Fig. 1 can be subdivided into three stages, i.e. the pre-stretching, the die’s uprising and the final stretching. During the pre-stretching, the metal skin is stretched horizontally to generate the pre-stress. Then the die rises and pushes the skin to generate a curved surface. Finally, the skin is pulled along the clamping direction to accomplish the final stretching. The simulation of the stretch-forming procedure can be solved by standard nonlinear finite element method. In this paper, we choose to use ABAQUS solver with central difference method and explicit
280
J.-H. Zhu et al. / Journal of Computational and Applied Mathematics 246 (2013) 278–288
Fig. 2. An illustrative definition and solution of the topology optimization problem.
time-integration, where the dynamic equation can be written as
¨ (t ) = M −1 · F (t ) − I (t ) u
(1)
where M is the mass matrix; F (t ) is nodal load vector at the moment t; I (t ) is the inertial force vector. The explicit timeintegration can be described as t + 12t
˙ u
t − 12t
= u˙
+
1t (t +1t ) + 1t (t ) 2
¨ (t ) u
t + 12t ˙ u(t +1t ) = u(t ) + 1t (t +1t ) u .
(2)
(3)
Based on the strain velocity ε˙ , the element strain increment dε is calculated. According to the material constitutive equation, the stress can be written as
σ (t +1t ) = f σ (t ) , dε .
(4)
Finally, the internal force I (t +1t ) is composed with the obtained stress. Replacing the current time t by t + 1t, the nonlinear problem is then solved iteratively. Using the obtained contact surface loads, Fig. 2 shows an illustrative definition and solution of the topology optimization problem. The working surface of the die is defined as the non-designable area, and the inner part is assigned as the design domain. To maximize the structural stiffness, the objective function can thus be expressed as a weighted function of the strain energies calculated from several typical load cases.
min :
Cw =
k l =1
w(l) C(l) =
k 1 l=1
2
T w(l) f(l) + G U(l)
(5)
where Cw is the weighted strain energy; C(l) is the strain energy of the l-th load case; w(l) is the weight factor of each load case; f(l) is the load vector and G is the self-weight load vector; U(l) is the corresponding nodal displacement vector; k is the total number of the load cases.
2.2. Topology optimization and casting constraints Moreover, to realize the casting constraints, a linear interpolation approach is applied in this paper. The design variables for the topology optimization are actually defined as the height of each column along the casting direction. It is further interpolated with the pseudo-densities of the related elements. Similar idea has been proposed by Gersborg and Andreasen [14], where the height of each column is interpolated with a Heaviside function. While in this paper, we propose to use an alternative linear interpolation to simplify the definition. The idea is illustrated in Fig. 3. Suppose the size of the element is 1, an artificial straight line is defined as the cutting line. The
J.-H. Zhu et al. / Journal of Computational and Applied Mathematics 246 (2013) 278–288
281
Fig. 3. The definition of the linear interpolation approach for the casting constraints.
elements in a column below the line are assigned as solid and the above as void. Only the elements cut by the line will be interpolated according to the divided dark areas. For a single column, we can have hi =
m
xi(j)
j =1
xi(j)
1; (⌈hi ⌉ − hi )2 ; H1 (hi ) = 1 − 4 1 h − ⌊hi ⌋ = H2 (hi ) = + i ; 4 2 (hi − ⌊hi ⌋)2 H3 (hi ) = ; 4 0;
j < ⌊hi ⌋ j = ⌊hi ⌋ j = ⌈hi ⌉
j = 1, 2, . . . , m
(6)
j = ⌈hi ⌉ + 1 j > ⌈hi ⌉ + 1
where hi indicates the height of the i-th column. ⌊hi ⌋ and ⌈hi ⌉ are the floor and ceiling functions of hi . H1 , H2 and H3 indicate 3 different dark areas. m is the total number of the elements in this column. xi(j) indicates the pseudo-density of the j-th element in the i-th column. The values of the pseudo-densities describe the solid-void material distribution in the topology optimization with a material interpolation model. Suppose E and Ke are the elastic modulus and the stiffness matrix of an element, E0 and Ke0 are the elastic modulus and the stiffness matrix of a solid element, x is the pseudo-density, we have E = P (x) E0
(7)
Ke = P (x) K e0
where P (x) is the so called penalty function forcing the pseudo-densities to vary towards 0 or 1 [2,3]. In this paper, the effect of penalty is useless as the pseudo-densities are assigned according to the height of the column. So we can simply define P (x) = x.
(8)
As a result, the topology optimization problem can be expressed as follows. The heights of the columns are defined as the design variables to describe the material distribution. The compliances, i.e. the global strain energies obtained from different load cases will be weighted and summed as the objective function. The total material cost is constrained to save the weight. find:
( h1 , h2 , . . . , hn )
min :
Cw =
k
w(l) C(l)
l =1
s.t.:
V ≤ V0
where V0 is the maximum allowed material cost.
(9)
282
J.-H. Zhu et al. / Journal of Computational and Applied Mathematics 246 (2013) 278–288
(a) The design domain, loads and boundary conditions.
(b) The optimized design without casing constraints.
(c) The optimized design with the linear interpolation.
Fig. 4. A simple example to show the effect of the casting constraint.
Fig. 5. The affected elements change from {1, 2, 3} to {2, 3, 4} when hi is an integer.
2.3. Sensitivity analysis According to Eqs. (6)–(8), the sensitivities of the strain energies with respect to hi can be calculated by the derivatives with respect to the pseudo-densities in the corresponding column. We can directly have
∂ C(l) ∂ C(l) ∂ H1 ∂ C(l) ∂ H2 ∂ C(l) ∂ H3 = + + . ∂ hi ∂ xi(⌊hi ⌋) ∂ hi ∂ xi(⌈hi ⌉) ∂ hi ∂ xi(⌈hi ⌉+1) ∂ hi
(10)
The derivative with respect to the pseudo-densities can be further derived as
T T ∂ U(l) ∂ f(l) + G U(l) + f(l) + G ∂ xi(j) ∂ xi(j) T T −1 ∂ K 1 ∂G = U(l) − f(l) + G K U(l) 2 ∂ xi(j) ∂ xi(j) 1 ∂ Ki(j) = GiT(j)0 Ui(j) − UiT(j) Ui(j) 2 ∂ xi(j) 1 T = Gi(j)0 Ui(j) − UiT(j) Ki(j)0 Ui(j)
1 ∂ C(l) = ∂ xi(j) 2
2
(11)
where Gi(j) and Ki(j) are the self-weight load vector and the stiffness matrix of the j-th element in the i-th column. Gi(j)0 and Ki(j)0 are the corresponding self-weight load vector and the stiffness matrix when the element is pure solid. A simple cantilever beam is designed to test the effect of the linear interpolation. The design domain, loads and boundary conditions are shown in Fig. 4(a). The problem is first solved without any casting constraints with the optimized design shown in Fig. 4(b). Later, the casting direction is assigned and the linear interpolation is implemented. The optimized design is shown in Fig. 4(c), where an obvious difference of the structural topology can be found. Since the linear interpolation function is a piecewise function, the continuities of the objective function and its derivatives have to be proved. The continuity of the function is easily maintained. To further inspect the continuity of the derivatives, we find the possible discontinuity only happens when hi is an integer. As the elements shown in Fig. 5, the previously affected elements shift here from {1, 2, 3} to {2, 3, 4}. As a result, the left derivative with respect to hi can be expressed as
∂ C(l) ∂ H1 ∂ C(l) ∂ H2 ∂ C(l) ∂ H3 ∂ C(l) = + + ∂ hi ∂ x1 ∂ hi x1 =1 ∂ x2 ∂ hi x2 =0.75 ∂ x3 ∂ hi x3 =0.25 =
∂ C(l) 1 ∂ C(l) 1 ∂ C(l) ·0+ · + · . ∂ x1 ∂ x2 2 ∂ x3 2
(12)
Note that in the left derivative, we have hi − ⌊hi ⌋ = 1
⌈hi ⌉ − hi = 0.
(13)
J.-H. Zhu et al. / Journal of Computational and Applied Mathematics 246 (2013) 278–288
283
Fig. 6. The elements selected from the optimized design and the change of the height.
Fig. 7. The variation of the global strain energy versus the height of the column.
The right derivative can be expressed similarly as
∂ C(l) ∂ H2 ∂ C(l) ∂ H3 ∂ C(l) ∂ C(l) ∂ H1 + + = ∂ hi ∂ x2 ∂ hi x2 =0.75 ∂ x3 ∂ hi x3 =0.25 ∂ x4 ∂ hi x4 =0 =
∂ C(l) 1 ∂ C(l) 1 ∂ C(l) · + · + · 0. ∂ x2 2 ∂ x3 2 ∂ x4
(14)
Note that in the right derivative, we have hi − ⌊ hi ⌋ = 0
⌈hi ⌉ − hi = 1.
(15)
In this way, we prove that the continuity of the derivative is strictly maintained. We can further show how the structural strain energy and the derivatives change over the height of a column in a numerical test. Several elements in the circle as shown in Fig. 6 are selected from the optimized design in Fig. 4(c). The height of the column rises by 3 elements. The procedure is subdivided into 300 short steps. The variation of the global strain energy is plotted in Fig. 7. The continuity and the differentiability of the global strain energy are well presented in this curve. As a result, we can conclude that the linear interpolation can equally ensure the continuity and convergence of the optimization. It also holds the advantage of simplicity compared with the complex and high nonlinear Heaviside function interpolation. The computational efficiency can also be significantly improved by using the height variables. On the one hand, the number of the height variables is much less than the number of pseudo-densities compared with the classical topology optimization. On the other hand, no extra design constraint is introduced to control the material distribution, which saves more computational costs during optimization solution. 3. Numerical simulation A certain type of skin stretch-forming system and the traditional design of the die are shown in Figs. 8 and 9. The maximum dimensions of the dual-curvature horizontal die are 2000 mm (Length) × 1000 mm (Width) × 480 mm (Height). The sizes of the metal skin are 2000 mm (Width) × 3000 mm (Length) × 2 mm (Thickness), which is clamped at the two ends. The die is discretized into fine hexahedron solid elements with the average size 20 mm. Meanwhile, the metal skin is discretized into quadrangular shell elements with the size 20 mm. The materials of the die and the aircraft skin are set to be carbon steel Q235A and 2024-O, respectively. The true stress-plastic strain curve of 2024-O is shown in Fig. 10. Coulomb Friction Model [15] is used to describe the friction between the die and skin, where the friction coefficient is assigned as 0.103. In this paper, the skin stretch-forming process consists of the so called pre-stretching and the die’s uprising. The skin is first clamped at the two ends and stretched horizontally until 80% of the yield limit. Then the die rises to push the skin. The metal sheet forms with the die gradually where the surface load varies during the forming process.
284
J.-H. Zhu et al. / Journal of Computational and Applied Mathematics 246 (2013) 278–288
Fig. 8. CAD model of the stretch-forming system.
Fig. 9. Traditional design of the stretch-forming die.
Fig. 10. 2024-O true stress-plastic strain curve.
Actually, the contact surface loads are design-dependent and relying on the stiffness of the die. In a strict optimization design procedure, the loads shall be updated with the design iterations. However, on account of the computational difficulties involved in the iterative nonlinear calculations, the updated surface loads as well as their design sensitivities are extremely computational expensive. As a result, we decide to use constant surface loads approximately. To validate the surface loads used in the topology optimization, we compared the surface loads obtained with three different stiffness configurations of the die, i.e. rigid, solid and the traditional design respectively. Several time steps indicated as C1 to C5 in Fig. 11 are selected as the key load cases where the nodal surface loads are local maximums on the curve. The maximum values and the distributions of the surface loads at different time steps are obtained and plotted in Figs. 11 and 12. In these Figures, we can find the distributions and the maximum values of the surface loads on the solid and the traditional dies are very similar with each other. Only the loads of the rigid die are different from those of the two elastic ones. Moreover, the stiffness of a reasonable optimized design has to lie somewhere between the solid and the traditional design. As a result, we decide to choose the 5 key load cases from the solid die as the approximate load cases for the topology optimization. 4. Topology optimization design of stretch-forming die In the initial model of the stretch-forming die, a working surface of 10 mm thickness is assigned as the non-designable domain. The surface loads and the self-weight are applied as multiple load cases with the bottom of the die clamped. 20% of the material cost allowed, the topology optimization is carried out with the weight ratio of five load cases ω1 :ω2 :ω3 :ω4 :ω5 = 1:1:1:2:2. The 4th load case is the extreme condition. Meanwhile the surface loads in the 5th load case are concentrated at the two sides of the working surface, which may lead to large deformation. We applied the linear interpolation with the casting direction assigned from the top to the bottom of the die. The optimized design is shown in Fig. 13.
J.-H. Zhu et al. / Journal of Computational and Applied Mathematics 246 (2013) 278–288
285
Fig. 11. Maximum nodal loads versus time steps in stretch-forming process.
Fig. 12. The distribution of the surface loads on the working surface when the die is (a) pure rigid; (b) pure solid; (c) the traditional design.
Later, we strengthen the working surface from 10 mm to 15 mm according to the engineering suggestions. Although the optimized design is much stiffer than the previous one, the material distribution is almost identical with Fig. 13. To inspect the effect of topology optimization, four different designs are analyzed and compared here. The maximum deformations and Von-Mises stresses are listed in Table 1. The four different designs are: (1) The traditional design as shown in Fig. 9. (2) The topological design as shown in Fig. 13. (3) Based on the optimized design in Fig. 13, the die is reconstructed in CAD as shown in Fig. 14.
286
J.-H. Zhu et al. / Journal of Computational and Applied Mathematics 246 (2013) 278–288
Fig. 13. Topological design of the die.
Table 1 The comparison of the four different designs. Designs
Volume fraction (%)
Maximum deformation
Von-Mises stress
(1) (2) (3) (4)
21.2 20 20.9 20.9
0.0755 mm 0.0482 mm 0.0519 mm 0.0334 mm
64.24 MPa N/A 23.80 MPa 20.85 MPa
Fig. 14. Directly reconstructed optimized result for practical engineering.
Fig. 15. Reconstructed optimized result with the thickness of the working surface strengthened to 15 mm.
(4) The working surface is strengthened to 15 mm. Some of the material is moved from the supporting structures to the surface to keep the weight identical with the previous one. The die is reconstructed in CAD as shown in Fig. 15. Compared with the traditional design (1), the topological design (2) can improve the stiffness of the die significantly. The stress is not provided here because it is not a pure solid-void design and the stress is not convincible. In the reconstructed design (3), the stiffness of the die is inevitably a little bit weaker than the optimized design (2). Compared with the traditional one, the maximum deformation and the maximum Von-Mises stress are obtained as 0.0519 mm and 23.80 MPa, which decrease by 31.3% and 63%, respectively. In the reconstructed design (4), the working surface is strengthened from 10 mm to 15 mm. The maximum deformation and the Von-Mises stress are 0.0334 mm and 20.85 MPa, which are 55.8% and 67.5% lower than the traditional design, respectively. The variation histories of the maximum deformation and stress level at different time steps are compared in Figs. 16 and 17. We can find the optimized designs have improved the stiffness and strength of the die significantly. Especially, in the topological design with the strengthened working surface of 15 mm is globally better than the traditional and the first optimized design.
J.-H. Zhu et al. / Journal of Computational and Applied Mathematics 246 (2013) 278–288
287
Fig. 16. Maximum deformations of different designs versus time steps.
Fig. 17. Maximum stress values of different designs versus time steps.
5. Conclusions The topology optimization of the aircraft skin stretch forming die is studied in this paper. Based on the stretch-forming simulation, the surface loads on the die are obtained and applied as the multiple load cases. Although the surface loads are design dependent and relying on the stiffness of the die, further analyses have proved that the loads can be approximately treated as constant in a certain range, which favors the iterative calculation of the optimization. As a result, we choose to use the surface loads of a solid die. Then the structural design of the die was carried out with the topology optimization method. We introduced a linear interpolation approach for the casting constraints. The material distribution is constrained on the prescribed casting direction. We further prove that the linear interpolation of the pseudo-densities can strictly maintain the continuity and the differentiability. Different optimized designs are analyzed and compared with the traditional design. The deformations and the stress levels of the optimized designs decrease significantly with the total weight maintained. Acknowledgments The first author’s work is supported by National Natural Science Foundation for Distinguished Young Scholar (10925212), National Natural Science Foundation of China (90916027, 11002113, 11172236), 973 Program (2011CB610304), the 111 Project (B07050) and NPU Foundation for Fundamental Research (NPU-FFR-JC20120229). References [1] M.P. Bendsoe, O. Sigmund, Topology Optimization: Theory, Methods, and Applications, Springer, Berlin, Heidelberg, New York, 2003. [2] M.P. Bendsoe, Optimal shape design as a material distribution problem, Structural Optimization 10 (1989) 193–202. [3] M. Zhou, G.I.N. Rozvany, The COC algorithm, part II: topological, geometry and generalized shape optimization, Computer Methods in Applied Mechanics and Engineering 89 (1991) 197–224.
288
J.-H. Zhu et al. / Journal of Computational and Applied Mathematics 246 (2013) 278–288
[4] J. París, F. Navarrina, I. Colominas, M. Casteleiro, Improvements in the treatment of stress constraints in structural topology optimization problems, Journal of Computational and Applied Mathematics 234 (2010) 2231–2238. [5] A. Erman Tekkaya, State-of-the-art of simulation of sheet metal forming, Journal of Materials Processing Technology 103 (2000) 14–22. [6] A.M. Yan, I. Klappka, Springback in stretch forming process of aeronautic panel production by finite element simulation, International Journal of Material Forming 1 (2008) 201–204. [7] S.H. Wang, Z.Y. Cai, M.Z. Li, FE simulation of shape accuracy using the Multi-Point Stretch-Forming process, International Journal of Materials and Product Technology 38 (2010) 223–236. [8] Z.Y. Jiang, S.Q. Li, J.M. Zeng, X.P. Liao, D.G. Yang, Numeral simulation of a new forming process: stretch forming combined with roll forming, Advanced Materials Research 189–193 (2011) 2399–2403. [9] J. Petersson, M. Patriksson, Topology optimization of sheets in contact by a subgradient method, International Journal for Numerical Methods In Engineering 40 (1997) 1295–1321. [10] N. Strömberg, A. Klarbring, Topology optimization of structures in unilateral contact, Structural and Multidisciplinary Optimization 41 (2010) 57–64. [11] J.J Sheu, C.H. Yang, A simplified column model for the automatic design of the stamping die structure, Journal of Materials Processing Technology 177 (2006) 109–113. [12] A. Nilsson, F. Birath, Topology optimization of a stamping die, in: The 9th International Conference on Numerical Methods in Industrial Forming Processes, Porto, Portugal, 2007. [13] D.K. Xu, Y.C. Tang, J. Chen, The topology optimization of stamping die of high-strength steel considered the bias load, Die and Mould Technology 4 (2009) 10–14. [14] A.R. Gersborg, C.S. Andreasen, An explicit parameterization for casting constraints in gradient driven topology optimization, Structure and Multidisciplinary Optimization 44 (2011) 875–881. [15] Z.Y. Cai, S.H. Wang, X.D. Xu, M.Z. Li, Numerical simulation for the multi-point stretch forming process of sheet metal, Journal of Materials Processing Technology 209 (2009) 396–407.
Contents lists available at SciVerse ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
Structural design of aircraft skin stretch-forming die using topology optimization Ji-Hong Zhu a,∗ , Xiao-Jun Gu a , Wei-Hong Zhang a , Pierre Beckers b a
Engineering Simulation and Aerospace Computing (ESAC), Northwestern Polytechnical University, Xi’an, Shaanxi 710072, China
b
LTAS-Infographie, University of Liège, Liège 4000, Belgium
article
info
Article history: Received 23 January 2012 Received in revised form 29 August 2012 Keywords: Stretch-forming die Light-weight design Topology optimization Linear interpolation Casting constraints
abstract This paper demonstrates a topology optimization method for the large scale stretchforming die design. The implementation of the topology optimization depends on the surface loads obtained from the numerical simulation of the sheet metal stretch-forming process. Typically, the design variables are defined based on a linear interpolation to generate directional structures satisfying the manufacturing requirement. To further validate the design procedure, mathematical derivations have proved that the design sensitivities are strictly continuous. More explanations and numerical tests are presented to show the variation of the surface loads versus the stiffness of the die. A stretch-forming die design example is solved on account of gravity and surface loads. The final solution is compared with the traditional design. The numerical results have shown that the topology design can improve the stiffness and strength of the stretchforming die significantly. © 2012 Elsevier B.V. All rights reserved.
1. Introduction In the past two decades, the techniques of topology optimization were developed to generate reasonable structural configurations by re-distributing the material in the design domains with the prescribed loads and boundary conditions [1–4]. Till now, one of the most challenging works involved in the research of topology optimization is the industrial application, where complex practical load cases have to be taken into account and different manufacturing constraints have to be satisfied. The skin stretch forming die discussed in this paper is one of the key equipments of aircraft manufacturing. It basically consists of a certain stretch working surface of specific form and the supporting structure that interconnects the working surface for its integrity. Nowadays, the development of large cargo and passenger aircraft has brought a series of new challenges. On the one hand, the size and weight of the stretch-forming dies increase significantly. As the traditional die is always designed extremely stiff, it results in excessive costs as well as difficulties in transportation and installation. On the other hand, the structural design of the die must ensure the accuracy of the stretch-forming process. As a result, the light-weight design of the large scaled stretch-forming die using the structural optimization methods becomes necessarily important. In this paper, the topology optimization method is used to design the supporting structures on account of the gravity of the die itself and the surface loads obtained from the numerical stretch forming procedure.
∗
Corresponding author. E-mail address: [email protected] (J.-H. Zhu).
0377-0427/$ – see front matter © 2012 Elsevier B.V. All rights reserved. doi:10.1016/j.cam.2012.09.001
J.-H. Zhu et al. / Journal of Computational and Applied Mathematics 246 (2013) 278–288
279
(b) The pre-stretching.
(c) The die’s uprising.
(a) A stretch-forming system and its working procedure (www.custompartnet.com).
(d) The final stretching.
Fig. 1. An illustration of stretch-forming procedure.
The structural optimization design of the stretch-forming die proposed in this paper is actually carried out as a twostep way. First, the surface loads on the die have to be numerically predicted with the performing of non-linear finite element analysis. Some typical and extreme load cases will be selected for the topology optimization. Currently, most of the existing works on the stretch-forming were concentrated on the numerical simulation of different kinds of forming process to improve the forming precision [5]. For example, Yan and Klappka [6] studied the spring-back behaviors of panel forming productions using the multi-point stretch forming technique. Wang et al. [7] carried out a series of numerical simulations. The shape error and non-uniformity in thickness distribution of the sheet part were analyzed. Jiang et al. [8] combined the stretch forming with the roll forming procedure. The non-linear Finite Element simulation had shown more benefits that better forming precision and lower stress level were obtained. Till now, the development of the numerical simulation of the metal sheet stretch forming procedure has been satisfied. As a result, we choose to use the explicit time-integration method in the platform of ABAQUS to obtain the different surface load cases. Second, based on the contact forces, topology optimization will be implemented to design the supporting structure satisfying the stiffness of the die and the requirement of weight-saving. Currently, the application of topology optimization to the contact problems is not a new issue. For example, Petersson and Patriksson [9] proposed a subgradient algorithm to solve the topology optimization of sheets in contact. The convergence and the computational efficiency of the optimization were in-depth discussed. Recently, in the work of Strömberg and Klarbring [10], the unilaterally constraints were simulated with contact conditions. Both 2D and 3D problems were optimized and reasonable results were obtained. In the existing works, the structural optimizations of the metal forming dies are mostly found designing the stamping die. For example, Sheu and Yang [11] predicted the pressures on the working surface of a stamping die. The inner structures were then designed by modifying the size and layout of the stiffeners. The optimized results turned out to be completely different from the traditional design, i.e. a uniform distribution of the stiffeners. Nilsson and Birath [12] simulated the variation of the surface loads during the stamping process by means of time integration, where the process of lifting and stamping are taken into account. Topology optimization was then applied to save the weight by maintaining the structure rigidity and strength. Similarly, Xu and Tang [13] developed the topology optimization procedure for the inner structure design of the stamping die based on the platform of LS-Dyna and Hyperworks, where 28% of the total weight is saved. Based on the previous works, a stretch-forming die of a large cargo aircraft is designed in this paper. The stretch-forming process is first simulated based on ABAQUS to obtain the surface loads which will be defined as numbers of weighted load cases. By deriving the design sensitivities and applying the constraints on the material distribution and the total material cost, topology optimization of the inner structure is then carried out to find the optimized structural layout. 2. Basic formulation 2.1. Stretch-forming simulation A typical stretch-forming procedure as illustrated in Fig. 1 can be subdivided into three stages, i.e. the pre-stretching, the die’s uprising and the final stretching. During the pre-stretching, the metal skin is stretched horizontally to generate the pre-stress. Then the die rises and pushes the skin to generate a curved surface. Finally, the skin is pulled along the clamping direction to accomplish the final stretching. The simulation of the stretch-forming procedure can be solved by standard nonlinear finite element method. In this paper, we choose to use ABAQUS solver with central difference method and explicit
280
J.-H. Zhu et al. / Journal of Computational and Applied Mathematics 246 (2013) 278–288
Fig. 2. An illustrative definition and solution of the topology optimization problem.
time-integration, where the dynamic equation can be written as
¨ (t ) = M −1 · F (t ) − I (t ) u
(1)
where M is the mass matrix; F (t ) is nodal load vector at the moment t; I (t ) is the inertial force vector. The explicit timeintegration can be described as t + 12t
˙ u
t − 12t
= u˙
+
1t (t +1t ) + 1t (t ) 2
¨ (t ) u
t + 12t ˙ u(t +1t ) = u(t ) + 1t (t +1t ) u .
(2)
(3)
Based on the strain velocity ε˙ , the element strain increment dε is calculated. According to the material constitutive equation, the stress can be written as
σ (t +1t ) = f σ (t ) , dε .
(4)
Finally, the internal force I (t +1t ) is composed with the obtained stress. Replacing the current time t by t + 1t, the nonlinear problem is then solved iteratively. Using the obtained contact surface loads, Fig. 2 shows an illustrative definition and solution of the topology optimization problem. The working surface of the die is defined as the non-designable area, and the inner part is assigned as the design domain. To maximize the structural stiffness, the objective function can thus be expressed as a weighted function of the strain energies calculated from several typical load cases.
min :
Cw =
k l =1
w(l) C(l) =
k 1 l=1
2
T w(l) f(l) + G U(l)
(5)
where Cw is the weighted strain energy; C(l) is the strain energy of the l-th load case; w(l) is the weight factor of each load case; f(l) is the load vector and G is the self-weight load vector; U(l) is the corresponding nodal displacement vector; k is the total number of the load cases.
2.2. Topology optimization and casting constraints Moreover, to realize the casting constraints, a linear interpolation approach is applied in this paper. The design variables for the topology optimization are actually defined as the height of each column along the casting direction. It is further interpolated with the pseudo-densities of the related elements. Similar idea has been proposed by Gersborg and Andreasen [14], where the height of each column is interpolated with a Heaviside function. While in this paper, we propose to use an alternative linear interpolation to simplify the definition. The idea is illustrated in Fig. 3. Suppose the size of the element is 1, an artificial straight line is defined as the cutting line. The
J.-H. Zhu et al. / Journal of Computational and Applied Mathematics 246 (2013) 278–288
281
Fig. 3. The definition of the linear interpolation approach for the casting constraints.
elements in a column below the line are assigned as solid and the above as void. Only the elements cut by the line will be interpolated according to the divided dark areas. For a single column, we can have hi =
m
xi(j)
j =1
xi(j)
1; (⌈hi ⌉ − hi )2 ; H1 (hi ) = 1 − 4 1 h − ⌊hi ⌋ = H2 (hi ) = + i ; 4 2 (hi − ⌊hi ⌋)2 H3 (hi ) = ; 4 0;
j < ⌊hi ⌋ j = ⌊hi ⌋ j = ⌈hi ⌉
j = 1, 2, . . . , m
(6)
j = ⌈hi ⌉ + 1 j > ⌈hi ⌉ + 1
where hi indicates the height of the i-th column. ⌊hi ⌋ and ⌈hi ⌉ are the floor and ceiling functions of hi . H1 , H2 and H3 indicate 3 different dark areas. m is the total number of the elements in this column. xi(j) indicates the pseudo-density of the j-th element in the i-th column. The values of the pseudo-densities describe the solid-void material distribution in the topology optimization with a material interpolation model. Suppose E and Ke are the elastic modulus and the stiffness matrix of an element, E0 and Ke0 are the elastic modulus and the stiffness matrix of a solid element, x is the pseudo-density, we have E = P (x) E0
(7)
Ke = P (x) K e0
where P (x) is the so called penalty function forcing the pseudo-densities to vary towards 0 or 1 [2,3]. In this paper, the effect of penalty is useless as the pseudo-densities are assigned according to the height of the column. So we can simply define P (x) = x.
(8)
As a result, the topology optimization problem can be expressed as follows. The heights of the columns are defined as the design variables to describe the material distribution. The compliances, i.e. the global strain energies obtained from different load cases will be weighted and summed as the objective function. The total material cost is constrained to save the weight. find:
( h1 , h2 , . . . , hn )
min :
Cw =
k
w(l) C(l)
l =1
s.t.:
V ≤ V0
where V0 is the maximum allowed material cost.
(9)
282
J.-H. Zhu et al. / Journal of Computational and Applied Mathematics 246 (2013) 278–288
(a) The design domain, loads and boundary conditions.
(b) The optimized design without casing constraints.
(c) The optimized design with the linear interpolation.
Fig. 4. A simple example to show the effect of the casting constraint.
Fig. 5. The affected elements change from {1, 2, 3} to {2, 3, 4} when hi is an integer.
2.3. Sensitivity analysis According to Eqs. (6)–(8), the sensitivities of the strain energies with respect to hi can be calculated by the derivatives with respect to the pseudo-densities in the corresponding column. We can directly have
∂ C(l) ∂ C(l) ∂ H1 ∂ C(l) ∂ H2 ∂ C(l) ∂ H3 = + + . ∂ hi ∂ xi(⌊hi ⌋) ∂ hi ∂ xi(⌈hi ⌉) ∂ hi ∂ xi(⌈hi ⌉+1) ∂ hi
(10)
The derivative with respect to the pseudo-densities can be further derived as
T T ∂ U(l) ∂ f(l) + G U(l) + f(l) + G ∂ xi(j) ∂ xi(j) T T −1 ∂ K 1 ∂G = U(l) − f(l) + G K U(l) 2 ∂ xi(j) ∂ xi(j) 1 ∂ Ki(j) = GiT(j)0 Ui(j) − UiT(j) Ui(j) 2 ∂ xi(j) 1 T = Gi(j)0 Ui(j) − UiT(j) Ki(j)0 Ui(j)
1 ∂ C(l) = ∂ xi(j) 2
2
(11)
where Gi(j) and Ki(j) are the self-weight load vector and the stiffness matrix of the j-th element in the i-th column. Gi(j)0 and Ki(j)0 are the corresponding self-weight load vector and the stiffness matrix when the element is pure solid. A simple cantilever beam is designed to test the effect of the linear interpolation. The design domain, loads and boundary conditions are shown in Fig. 4(a). The problem is first solved without any casting constraints with the optimized design shown in Fig. 4(b). Later, the casting direction is assigned and the linear interpolation is implemented. The optimized design is shown in Fig. 4(c), where an obvious difference of the structural topology can be found. Since the linear interpolation function is a piecewise function, the continuities of the objective function and its derivatives have to be proved. The continuity of the function is easily maintained. To further inspect the continuity of the derivatives, we find the possible discontinuity only happens when hi is an integer. As the elements shown in Fig. 5, the previously affected elements shift here from {1, 2, 3} to {2, 3, 4}. As a result, the left derivative with respect to hi can be expressed as
∂ C(l) ∂ H1 ∂ C(l) ∂ H2 ∂ C(l) ∂ H3 ∂ C(l) = + + ∂ hi ∂ x1 ∂ hi x1 =1 ∂ x2 ∂ hi x2 =0.75 ∂ x3 ∂ hi x3 =0.25 =
∂ C(l) 1 ∂ C(l) 1 ∂ C(l) ·0+ · + · . ∂ x1 ∂ x2 2 ∂ x3 2
(12)
Note that in the left derivative, we have hi − ⌊hi ⌋ = 1
⌈hi ⌉ − hi = 0.
(13)
J.-H. Zhu et al. / Journal of Computational and Applied Mathematics 246 (2013) 278–288
283
Fig. 6. The elements selected from the optimized design and the change of the height.
Fig. 7. The variation of the global strain energy versus the height of the column.
The right derivative can be expressed similarly as
∂ C(l) ∂ H2 ∂ C(l) ∂ H3 ∂ C(l) ∂ C(l) ∂ H1 + + = ∂ hi ∂ x2 ∂ hi x2 =0.75 ∂ x3 ∂ hi x3 =0.25 ∂ x4 ∂ hi x4 =0 =
∂ C(l) 1 ∂ C(l) 1 ∂ C(l) · + · + · 0. ∂ x2 2 ∂ x3 2 ∂ x4
(14)
Note that in the right derivative, we have hi − ⌊ hi ⌋ = 0
⌈hi ⌉ − hi = 1.
(15)
In this way, we prove that the continuity of the derivative is strictly maintained. We can further show how the structural strain energy and the derivatives change over the height of a column in a numerical test. Several elements in the circle as shown in Fig. 6 are selected from the optimized design in Fig. 4(c). The height of the column rises by 3 elements. The procedure is subdivided into 300 short steps. The variation of the global strain energy is plotted in Fig. 7. The continuity and the differentiability of the global strain energy are well presented in this curve. As a result, we can conclude that the linear interpolation can equally ensure the continuity and convergence of the optimization. It also holds the advantage of simplicity compared with the complex and high nonlinear Heaviside function interpolation. The computational efficiency can also be significantly improved by using the height variables. On the one hand, the number of the height variables is much less than the number of pseudo-densities compared with the classical topology optimization. On the other hand, no extra design constraint is introduced to control the material distribution, which saves more computational costs during optimization solution. 3. Numerical simulation A certain type of skin stretch-forming system and the traditional design of the die are shown in Figs. 8 and 9. The maximum dimensions of the dual-curvature horizontal die are 2000 mm (Length) × 1000 mm (Width) × 480 mm (Height). The sizes of the metal skin are 2000 mm (Width) × 3000 mm (Length) × 2 mm (Thickness), which is clamped at the two ends. The die is discretized into fine hexahedron solid elements with the average size 20 mm. Meanwhile, the metal skin is discretized into quadrangular shell elements with the size 20 mm. The materials of the die and the aircraft skin are set to be carbon steel Q235A and 2024-O, respectively. The true stress-plastic strain curve of 2024-O is shown in Fig. 10. Coulomb Friction Model [15] is used to describe the friction between the die and skin, where the friction coefficient is assigned as 0.103. In this paper, the skin stretch-forming process consists of the so called pre-stretching and the die’s uprising. The skin is first clamped at the two ends and stretched horizontally until 80% of the yield limit. Then the die rises to push the skin. The metal sheet forms with the die gradually where the surface load varies during the forming process.
284
J.-H. Zhu et al. / Journal of Computational and Applied Mathematics 246 (2013) 278–288
Fig. 8. CAD model of the stretch-forming system.
Fig. 9. Traditional design of the stretch-forming die.
Fig. 10. 2024-O true stress-plastic strain curve.
Actually, the contact surface loads are design-dependent and relying on the stiffness of the die. In a strict optimization design procedure, the loads shall be updated with the design iterations. However, on account of the computational difficulties involved in the iterative nonlinear calculations, the updated surface loads as well as their design sensitivities are extremely computational expensive. As a result, we decide to use constant surface loads approximately. To validate the surface loads used in the topology optimization, we compared the surface loads obtained with three different stiffness configurations of the die, i.e. rigid, solid and the traditional design respectively. Several time steps indicated as C1 to C5 in Fig. 11 are selected as the key load cases where the nodal surface loads are local maximums on the curve. The maximum values and the distributions of the surface loads at different time steps are obtained and plotted in Figs. 11 and 12. In these Figures, we can find the distributions and the maximum values of the surface loads on the solid and the traditional dies are very similar with each other. Only the loads of the rigid die are different from those of the two elastic ones. Moreover, the stiffness of a reasonable optimized design has to lie somewhere between the solid and the traditional design. As a result, we decide to choose the 5 key load cases from the solid die as the approximate load cases for the topology optimization. 4. Topology optimization design of stretch-forming die In the initial model of the stretch-forming die, a working surface of 10 mm thickness is assigned as the non-designable domain. The surface loads and the self-weight are applied as multiple load cases with the bottom of the die clamped. 20% of the material cost allowed, the topology optimization is carried out with the weight ratio of five load cases ω1 :ω2 :ω3 :ω4 :ω5 = 1:1:1:2:2. The 4th load case is the extreme condition. Meanwhile the surface loads in the 5th load case are concentrated at the two sides of the working surface, which may lead to large deformation. We applied the linear interpolation with the casting direction assigned from the top to the bottom of the die. The optimized design is shown in Fig. 13.
J.-H. Zhu et al. / Journal of Computational and Applied Mathematics 246 (2013) 278–288
285
Fig. 11. Maximum nodal loads versus time steps in stretch-forming process.
Fig. 12. The distribution of the surface loads on the working surface when the die is (a) pure rigid; (b) pure solid; (c) the traditional design.
Later, we strengthen the working surface from 10 mm to 15 mm according to the engineering suggestions. Although the optimized design is much stiffer than the previous one, the material distribution is almost identical with Fig. 13. To inspect the effect of topology optimization, four different designs are analyzed and compared here. The maximum deformations and Von-Mises stresses are listed in Table 1. The four different designs are: (1) The traditional design as shown in Fig. 9. (2) The topological design as shown in Fig. 13. (3) Based on the optimized design in Fig. 13, the die is reconstructed in CAD as shown in Fig. 14.
286
J.-H. Zhu et al. / Journal of Computational and Applied Mathematics 246 (2013) 278–288
Fig. 13. Topological design of the die.
Table 1 The comparison of the four different designs. Designs
Volume fraction (%)
Maximum deformation
Von-Mises stress
(1) (2) (3) (4)
21.2 20 20.9 20.9
0.0755 mm 0.0482 mm 0.0519 mm 0.0334 mm
64.24 MPa N/A 23.80 MPa 20.85 MPa
Fig. 14. Directly reconstructed optimized result for practical engineering.
Fig. 15. Reconstructed optimized result with the thickness of the working surface strengthened to 15 mm.
(4) The working surface is strengthened to 15 mm. Some of the material is moved from the supporting structures to the surface to keep the weight identical with the previous one. The die is reconstructed in CAD as shown in Fig. 15. Compared with the traditional design (1), the topological design (2) can improve the stiffness of the die significantly. The stress is not provided here because it is not a pure solid-void design and the stress is not convincible. In the reconstructed design (3), the stiffness of the die is inevitably a little bit weaker than the optimized design (2). Compared with the traditional one, the maximum deformation and the maximum Von-Mises stress are obtained as 0.0519 mm and 23.80 MPa, which decrease by 31.3% and 63%, respectively. In the reconstructed design (4), the working surface is strengthened from 10 mm to 15 mm. The maximum deformation and the Von-Mises stress are 0.0334 mm and 20.85 MPa, which are 55.8% and 67.5% lower than the traditional design, respectively. The variation histories of the maximum deformation and stress level at different time steps are compared in Figs. 16 and 17. We can find the optimized designs have improved the stiffness and strength of the die significantly. Especially, in the topological design with the strengthened working surface of 15 mm is globally better than the traditional and the first optimized design.
J.-H. Zhu et al. / Journal of Computational and Applied Mathematics 246 (2013) 278–288
287
Fig. 16. Maximum deformations of different designs versus time steps.
Fig. 17. Maximum stress values of different designs versus time steps.
5. Conclusions The topology optimization of the aircraft skin stretch forming die is studied in this paper. Based on the stretch-forming simulation, the surface loads on the die are obtained and applied as the multiple load cases. Although the surface loads are design dependent and relying on the stiffness of the die, further analyses have proved that the loads can be approximately treated as constant in a certain range, which favors the iterative calculation of the optimization. As a result, we choose to use the surface loads of a solid die. Then the structural design of the die was carried out with the topology optimization method. We introduced a linear interpolation approach for the casting constraints. The material distribution is constrained on the prescribed casting direction. We further prove that the linear interpolation of the pseudo-densities can strictly maintain the continuity and the differentiability. Different optimized designs are analyzed and compared with the traditional design. The deformations and the stress levels of the optimized designs decrease significantly with the total weight maintained. Acknowledgments The first author’s work is supported by National Natural Science Foundation for Distinguished Young Scholar (10925212), National Natural Science Foundation of China (90916027, 11002113, 11172236), 973 Program (2011CB610304), the 111 Project (B07050) and NPU Foundation for Fundamental Research (NPU-FFR-JC20120229). References [1] M.P. Bendsoe, O. Sigmund, Topology Optimization: Theory, Methods, and Applications, Springer, Berlin, Heidelberg, New York, 2003. [2] M.P. Bendsoe, Optimal shape design as a material distribution problem, Structural Optimization 10 (1989) 193–202. [3] M. Zhou, G.I.N. Rozvany, The COC algorithm, part II: topological, geometry and generalized shape optimization, Computer Methods in Applied Mechanics and Engineering 89 (1991) 197–224.
288
J.-H. Zhu et al. / Journal of Computational and Applied Mathematics 246 (2013) 278–288
[4] J. París, F. Navarrina, I. Colominas, M. Casteleiro, Improvements in the treatment of stress constraints in structural topology optimization problems, Journal of Computational and Applied Mathematics 234 (2010) 2231–2238. [5] A. Erman Tekkaya, State-of-the-art of simulation of sheet metal forming, Journal of Materials Processing Technology 103 (2000) 14–22. [6] A.M. Yan, I. Klappka, Springback in stretch forming process of aeronautic panel production by finite element simulation, International Journal of Material Forming 1 (2008) 201–204. [7] S.H. Wang, Z.Y. Cai, M.Z. Li, FE simulation of shape accuracy using the Multi-Point Stretch-Forming process, International Journal of Materials and Product Technology 38 (2010) 223–236. [8] Z.Y. Jiang, S.Q. Li, J.M. Zeng, X.P. Liao, D.G. Yang, Numeral simulation of a new forming process: stretch forming combined with roll forming, Advanced Materials Research 189–193 (2011) 2399–2403. [9] J. Petersson, M. Patriksson, Topology optimization of sheets in contact by a subgradient method, International Journal for Numerical Methods In Engineering 40 (1997) 1295–1321. [10] N. Strömberg, A. Klarbring, Topology optimization of structures in unilateral contact, Structural and Multidisciplinary Optimization 41 (2010) 57–64. [11] J.J Sheu, C.H. Yang, A simplified column model for the automatic design of the stamping die structure, Journal of Materials Processing Technology 177 (2006) 109–113. [12] A. Nilsson, F. Birath, Topology optimization of a stamping die, in: The 9th International Conference on Numerical Methods in Industrial Forming Processes, Porto, Portugal, 2007. [13] D.K. Xu, Y.C. Tang, J. Chen, The topology optimization of stamping die of high-strength steel considered the bias load, Die and Mould Technology 4 (2009) 10–14. [14] A.R. Gersborg, C.S. Andreasen, An explicit parameterization for casting constraints in gradient driven topology optimization, Structure and Multidisciplinary Optimization 44 (2011) 875–881. [15] Z.Y. Cai, S.H. Wang, X.D. Xu, M.Z. Li, Numerical simulation for the multi-point stretch forming process of sheet metal, Journal of Materials Processing Technology 209 (2009) 396–407.