A Method for Optimal Blank Shape Determination in Sheet Metal Forming Based on Numerical Simulations

Nowadays most production processes are designed to increase economic efficiency while, at the same time, not influencing product quality, especially when dealing with large serial production. In order to improve the forming process from the economic point of view, the main objective is to minimize the number of operation phases involved in the forming process. The cutting phase, as the final operation, is the phase most often targeted for elimination. Along with less material required for product production, the elimination of the cutting phase may also reduce the occurrence of manufacturing defects during the forming process, such as tearing and wrinkling. In some cases, however, the sheet cutting phase is unavoidable since excess material under the blankholder is required to achieve proper holding of the sheet metal [1]. When the cutting phase is eliminated, the accuracy of the product edge geometry should be provided by proper blank shape geometry. In most cases, the proper blank shape is determined experimentally by a trial and error procedure which requires several try-outs, causing the forming process design to be extremely time consuming and costly. Nowadays, this can be done virtually, based on computer simulations of the technological process under consideration, see, for example, [1] to [4]. Besides choosing a proper numerical approach and computational technique [5] and [6], advanced constitutive modelling [2], [4], [7] and [8] and proper material characterisation [4] and [7] to [9] is crucial for the computer simulation to be physically objective and trustful. Although, in contrast to direct analysis, a Finite Element (FE) inverse analysis approach, see [10] to [14], could possibly be used in sheet metal forming simulation in order to reduce the computer time consumption, such an approach is not recommended because it results in less accurate strain–stress state determination. This, in turn, can considerably influence the subsequent springback analysis, which is a key element in tool design optimisation. This paper describes a blank shape optimisation method, which allows determination of a blank shape such that the edge geometry deviation of the produced product with respect to the specified product geometry is obtained within a given tolerance. In this method, the optimal blank shape is determined in an iterative way, starting with an initial blank of approximately adequate shape, which over a series of subsequent iterations gets automatically adjusted to meet the tolerance criterion. Following the principal idea, which is a gradual changing of the given blank shape based on the computed product geometry and manifested edge deviation, the product geometry must be determined in each iteration by a corresponding forming process simulation, considering the blank shape as determined in the previous iteration. The method is designed in such a way that it enables optimal blank shape determination of products having a general 3D shape. That the method is capable of tackling rather complex product geometries efficiently is demonstrated in Section 4, where results of the numerical optimisation process are also experimentally validated. A Method for Optimal Blank Shape Determination in Sheet Metal Forming Based on Numerical Simulations Mole, N. – Cafuta, G. – Štok, B. Nikolaj Mole1 – Gašper Cafuta2 – Boris Štok1,* 1 University of Ljubljana, Faculty of Mechanical Engineering, Slovenia 2 Cimos d.d., Automotive industry, Slovenia


INTRODUCTION
Nowadays most production processes are designed to increase economic efficiency while, at the same time, not influencing product quality, especially when dealing with large serial production.In order to improve the forming process from the economic point of view, the main objective is to minimize the number of operation phases involved in the forming process.The cutting phase, as the final operation, is the phase most often targeted for elimination.Along with less material required for product production, the elimination of the cutting phase may also reduce the occurrence of manufacturing defects during the forming process, such as tearing and wrinkling.In some cases, however, the sheet cutting phase is unavoidable since excess material under the blankholder is required to achieve proper holding of the sheet metal [1].
When the cutting phase is eliminated, the accuracy of the product edge geometry should be provided by proper blank shape geometry.In most cases, the proper blank shape is determined experimentally by a trial and error procedure which requires several try-outs, causing the forming process design to be extremely time consuming and costly.Nowadays, this can be done virtually, based on computer simulations of the technological process under consideration, see, for example, [1] to [4].Besides choosing a proper numerical approach and computational technique [5] and [6], advanced constitutive modelling [2], [4], [7] and [8] and proper material characterisation [4] and [7] to [9] is crucial for the computer simulation to be physically objective and trustful.Although, in contrast to direct analysis, a Finite Element (FE) inverse analysis approach, see [10] to [14], could possibly be used in sheet metal forming simulation in order to reduce the computer time consumption, such an approach is not recommended because it results in less accurate strain-stress state determination.This, in turn, can considerably influence the subsequent springback analysis, which is a key element in tool design optimisation.
This paper describes a blank shape optimisation method, which allows determination of a blank shape such that the edge geometry deviation of the produced product with respect to the specified product geometry is obtained within a given tolerance.In this method, the optimal blank shape is determined in an iterative way, starting with an initial blank of approximately adequate shape, which over a series of subsequent iterations gets automatically adjusted to meet the tolerance criterion.Following the principal idea, which is a gradual changing of the given blank shape based on the computed product geometry and manifested edge deviation, the product geometry must be determined in each iteration by a corresponding forming process simulation, considering the blank shape as determined in the previous iteration.
The method is designed in such a way that it enables optimal blank shape determination of products having a general 3D shape.That the method is capable of tackling rather complex product geometries efficiently is demonstrated in Section 4, where results of the numerical optimisation process are also experimentally validated.

REVIEW OF SOME APPROACHES TO BLANK SHAPE OPTIMISATION
The problem of finding an appropriate blank shape, which would allow the production of a formed product with the edge geometry meeting the geometry tolerance requirements, as specified by the design, has been addressed by many authors.In the section, we give a brief review of different numerical approaches developed in this regard.
An optimisation method, where the initial blank shape approximation is determined by inverse FEM simulation and the blank shape is optimised iteratively afterwards using direct FEM computer simulations, is presented in [15] and [16] by Azaouzi et al.The applied blank shape modification consists of displacing, in each iteration, the blank edge nodes in the direction opposite to the geometry deviation obtained by comparing the computed and reference product geometry.After the new blank shape is determined, its domain is remeshed automatically using triangular FEs.Naceur et al. [17] presented a blank shape optimisation approach that is based on the coupling between the inverse approach used for the forming simulation and an evolutionary algorithm.Their goal was to minimize the size of the blank shape and still ensure that the product is made without tearing the sheet metal.Park et al. [18] introduced a deformation path method.The Ideal Forming direct inverse design method, see [13] and [14], was used to determine initial blank shape, which was further improved by an iterative procedure.In [19], Yeh et al. use the inverse true strain method (TSM) to obtain an initial approximation of the blank shape, whereas for the optimal blank shape determination a method based on adaptive-network-based fuzzy inference system (ANFIS) is applied.To achieve satisfactory results on the optimised initial blank shape, at least for the case considered in the paper, 200 learning cycles, each requiring a complete computer simulation of the corresponding forming process, had to be used in building the appropriate ANFIS database.Another approach, where the abductive network is used to predict the optimal blank shape for forming an elliptic cup using FEM computer simulations, is presented by Lin and Kwan in [20].In the optimisation procedure, the initially elliptic blank shape, with its geometry being specified in a polar coordinate system by four characteristic points, is subject to change by considering adequate variation of the respective points' radial coordinates, while keeping their angular coordinates fixed.An iterative method, where the product geometry is also computed using inverse FEM while the blank shape correction is made in each iteration manually, is introduced by Parsa and Pournia in [21].Due to the application of the inverse analysis approach, the computational time is significantly reduced, but it is achieved at the expense of loss of accuracy in the product geometry determination.
The objective of the blank shape optimisation presented in [22] by Pegada et al. is to determine the blank shape that allows forming of a circular cup of uniform height, where the respective material properties are assumed to be orthotropic.In each iteration, considering the established deviation in the cup height, the blank shape is adapted proportionally, assuming a value of 0.75 as an adequate scaling factor to obtain satisfactory convergence to the method.In the method introduced by Son and Shim in [23], the blank shape correction is made in the direction opposite to the geometry deviation obtained by comparing the computed and reference product geometry.Furthermore, correction of each edge point is computed by applying a scaling factor between 0.5 and 0.9, its actual value being defined by a ratio of the initial velocity to the total deformation path length.In [24], Hamammi et al. propose a method where the blank shape correction is made in the direction described by the positions of each node lying on the blank edge at the beginning and at the end of the forming process.Another method, where the correction of the blank shape is also based on taking the displacement path of the product edge nodes into account, is proposed by Fazli and Arezoo in [25].They proved that their method is slightly better in terms of accuracy and also in convergence than the previous two methods, [23] and [24].
The optimal blank shape can also be obtained by parameterisation of the blank geometry and using a sequential programming method for finding the optimal set of parameters, which is elaborated in [26] by Sattari et al.Similarly, in [27], Padmanabhana et al. investigate blank shape optimisation using parameterisation of the blank shape carried out by using parametric NURBS curves and the blank shape correction based on the control points displacement.The convergence of the latter method can be further improved by including the sensitivity analysis shown in [28] by Shim et al.
In principle, all those approaches to blank shape optimisation have in common a sequence of rather similar logical steps, which can be roughly summarised by the flow chart in Fig. 3.The approach we are going to propose in this paper does not differ in this regard.
The iterative approach we present here is implemented in such a way that the user must provide an approximately determined blank shape, which is taken as the initial blank shape for the optimisation procedure.As demonstrated by the study case, elaborated in Section 4, the proposed numerical approach is efficient enough that it does not require a very accurate determination of the initial blank shape, which is the case in [15], [16], [18] and [19].By using different strategies the user can still provide a better approximation for the initial blank shape, which is certainly advantageous.One thing that is common to all the methods presented in the review is that, in each iteration, a remeshing of the blank is required after blank shape correction.In our case, on the other hand, the mesh element topology is kept unchanged through the entire iterative procedure, while the corresponding FE mesh nodal points' adjustment in the iteration is performed in accordance with the displacement field obtained by solving an associate linear elastic boundary value problem.In this boundary value problem, the blank with the actual shape geometry before correction is subjected to prescribed boundary displacements, the imposed boundary displacements taken equal the required blank edge geometry correction in the iteration.With the blank FE mesh correction dealt with in this way, any type of FEs can be used and no sophisticated remeshing methods are needed.

MATHEMATICAL, MODELLING AND PHYSICAL PRELIMINARIES
In order to provide a corresponding mathematical framework for the numerical analysis, in this section we give some definitions and conventions on the geometry terminology used, which is followed by a determination of some related point topology quantities, such as surface and edge normal vectors' determination, and, finally, an approximation for analytical surface reconstruction.Some modelling and simulation assumptions are given in the last part of the section in order to focus our investigation on the key elements of the proposed numerical procedure for the blank shape optimisation described in Section 3.

Geometry and Point Topology Terminology
In this paper, we adopt the same geometry terminology as introduced by Cafuta et al. in [29].The product geometry prescribed by the design engineer is referred as the "reference geometry", whereas the product geometry computed in the simulation will be referred as the "actual geometry".Considering that the forming process simulations are performed on the basis of FEM, we are actually dealing with discretised geometries.In accordance with the notation G , introduced for a general surface point topology definition, we will use notations G ref and G act to specify point topologies related to the reference and actual product geometry, respectively.Similarly, notation G bl will be used to denote point topology related to the sheet blank geometry.Furthermore, Γ ref , Γ act and Γ bl will be used to specify the associated edge point topologies notations.All those surface and edge topologies will exclusively refer to the sheet's mid-surface.In addition, to support a numerical procedure of automatic sheet blank geometry adjustment, an auxiliary virtual surface with its point topology notation being Σ G will be generated, emanating from the reference product edge Γ ref considering its surface topological properties.The surface Σ G having properties Γ ref ⊂ Σ G may be considered as a target surface which is to be attained iteratively by the edge points Γ P i k of the actual product Γ act k as closely as possible.Since, in the numerical procedure, the unit normal vectors n i at point P i perpendicular either to the surface or to the edge will also be referenced, it is reasonable to make a symbolic distinction between them.Accordingly, the notations n i and Σ n i will be used for the surface normal vectors, whereas Γ n i will be used for the edge normal vectors.Similarly, to make a distinction between the points appertaining either to the product surfaces, G ref and G act , or to the auxiliary surface Σ G , the notations P i and Σ P i will be used with respective position vectors being P i and Σ P i .If point P i is located on the edge Γ, it will be denoted as Γ P i and its position specified by vector Γ P i .

Surface Normal Vector Determination
The accuracy of the surface normal vector determination at points of discretely given surface geometries plays an important, if not crucial, role in attaining convergence and speeding up the convergence rate of the entire optimisation procedure.A general surface point topology G given, we apply in this paper a two-step procedure in which for a point of interest, say P i , two auxiliary vectors, ′ n i and ′′ n i , are computed, respectively.
In the first step, the auxiliary vector ′ n i at point P i is determined, considering FE surface discretisation.
The following equation is used: In the second step, the normal vector n i is determined considering the analytically defined smooth surface S i through point P i .This surface is obtained by a corresponding interpolation through a point set Π i , built from points P j i j ( =1, 2,...,9) of the FE mesh.Those points are chosen from a sub-domain area surrounding the considered point P P i i ≡ 1 on a closest point's basis, which is applicable regardless of whether point P P i i ≡ 1 belongs to the interior of the FE mesh or to its boundary.

Mathematical representation of the surface S i , F x y z
i ( , , ) = 0 , is built by considering coordinates of the respective points of the set Π i .To avoid round-off error due to possible computing with large numbers, it is advantageous to have the surface S i representation defined with respect to a particular local coordinate system, say ( , , )    x y z , having its origin at point P P where the base vector e z = (0,0,1) defines the global coordinate axis z.
The rotation of the global coordinate system (x, y, z) to the local coordinate system ( , , )    x y z is given by the rotational transformation matrix ℜ( , ) r θ :    x y z , the coordinates of the position vector P i j are mapped to the local coordinate system in accordance with: With the points set  Π i consisting of nine points  P i j (Fig. 1b), it is appropriate to approximate surface S i in accordance with: where nine coefficients a mn are determined such that the surface S i is interpolated through the points  P i j .
Finally, the auxiliary vector ′′ n i normal to the surface f x y z i ( , , ) = 0    at point  P i 1 (0,0,0) can be determined by the gradient operator: By applying to the components of the above vector, given in the local coordinate system, the inverse mapping ℜ −1 ( , ) r θ , the components of the normal vector n i to the surface S i at point P i in the global coordinate system are obtained such that: If the point, at which the normal vector is to be computed, is positioned on the symmetry plane, the normal vector n i is computed by the same procedure considering the mirror elements over the symmetry plane.
The abovedescribed procedure of the surface normal vectors' determination is general.In Section 3, it will be carried out with reference to the reference and actual product geometry, G ref and G act , as well as to the auxiliary surface Σ G .

Edge Normal Vector Determination
We will refer generally to a surface point topology G and associate edge point topology Γ also in the determination of the unit vectors normal to the edge of a bounded surface.At a point of interest on the boundary Γ, say Γ P i , the edge normal vector Γ n i will be determined considering respective surface and boundary curvatures.While the surface curvature is characterised by the respective surface normal vector n i , with its determination being described in Section 2.2, the boundary curvature can be characterised by a unit vector Γ t i tangential to the surface boundary at point Γ P i .
The tangential vector Γ t i can be correspondingly approximated considering actual discretisation of the boundary Γ in the closest vicinity of the point Γ P i .With points Γ P i−1 and Γ P i+1 being the points adjacent to point Γ P i (see Fig. 2) the following equation may be applied: Finally, by evaluating the vector product: the edge normal vector Γ n i , defined in the global coordinate system, can be determined.The corresponding graphical representation, with vector Γ n i lying in the tangential plane to surface S i and pointing to its exterior, is given in Fig. 2. If the point, at which the edge normal vector is to be computed, is positioned on the symmetry plane, in order to achieve geometric symmetry conditions, the tangential vector Γ t i is defined by a unit vector normal to the symmetry plane.The edge normal vector Γ n i is afterwards computed by the same procedure.

Analytical Surface Reconstruction
Assuming that the considered FE surface discretisations are all based on a quadrilateral surface element meshing (see Fig. 1), we will develop, in this section, an approximation to the analytical surface reconstruction of a single quadrilateral surface element, considering its particular topological properties.
Let a quadrilateral surface S i be discretely defined by its nodal points P j i j ( =1, 2,3, 4) and respective surface normal vectors n i j at those points, determined following the procedure described in Section 2.2.This set of data represents twelve independent parameters, which can be used in analytical surface reconstruction of surface S i , F x y z i ( , , ) = 0 .This can be achieved by considering the following functional form: 3 3 ) = 0

F x y z z a a x a y a x y a x a y a x y a
with twelve coefficients a m m i ( =1, 2,...,12) to be determined.By imposing that the above surface representation meets given requirements at four interpolation points P i j the following system of linear equations is obtained: whose solution yields the value of the coefficients a m i .Here, it should be emphasised that by fulfilling the nodal surface properties in the above way element after element C 1 continuity at discrete points of the surface S i is achieved, which is significant for the quality of the overall surface approximation.

Modelling and Simulation Assumptions
When the blank shape optimisation is based on the results of a corresponding computer simulation of the considered forming process, it is exceedingly important that the springback numerical estimation should be as reliable as possible.In this regard, the resulting stress-strain state in the formed product prior to the springback, as well as the established change in the sheet metal thickness itself, are crucial for determining the extent of the manifested springback.To attain reliability of the computed results, the process parameters, such as the sheet material behaviour, tools' kinematics, sheet blank shape geometry, the clearance between the punch and the die, the blank-holder force and tribological conditions between the surfaces in contact, should be considered as realistically as possible in the simulation.For the purpose of this investigation, let us assume that all above parameters, except sheet blank shape geometry, are defined adequately and after the computer simulation yield a product of a shape, whose deviation from the shape of the reference product, when measured in the direction perpendicular to the product surface, is within the prescribed tolerance (see Cafuta et al. [29]).A possible deviation in the circumference contour, larger than is allowed, is subject to a corresponding adjustment in the sheet blank shape geometry, which is actually the topic of this paper.
Accordingly, we also assume that all issues associated with FEM, such as adequate choice of the finite element type and the integration of the constitutive equations, are not disputable.

NUMERICAL IMPLEMENTATION OF THE PROPOSED BLANK SHAPE OPTIMISATION
The proposed blank shape optimisation method consists basically of a sequence of steps executed iteratively following the flow chart in Fig. 3. Along with the initial blank shape geometry G bl 0 and reference product geometry G ref provided, a complete FE model with the tool geometry and forming process parameters (blank-holder forces, friction conditions, tool movement, etc.) must be specified in order to start the procedure.Since the initial blank shape geometry G bl 0 plays an important role in attaining convergence and computational efficiency of the described blank shape optimisation procedure, this geometry should be determined in a way that the edge geometry Γ act 0 of the product, obtained by the corresponding forming process simulation under given process conditions, does not experience too excessive a deviation with respect to the reference edge geometry Γ ref .Before starting the procedure the auxiliary surface Σ G is determined as explained in Section 3.1.

Fig. 3. Flow chart of the iterative procedure
The iterative procedure essentially consists of performing in iteration, say k th iteration, three consecutive steps.In the first step, the complete technological process, consisting of the forming stage with the actual blank shape geometry G bl k considered and subsequent springback stage occurring after removal of the tools, is simulated.This is followed by computing the deviation between the actual product edge geometry Γ act k and its reference edge geometry Γ ref , measured with respect to the auxiliary surface Σ G .The edge deviation computation is performed in each iteration following the procedure described in Section 3.2.In the last step, the actual blank shape geometry G bl k is either confirmed as appropriate or, if the considered edge deviation being found greater than the prescribed tolerance, its further correction and adjustment to G bl k +1 is required (see Section 3.3).Fulfilment of the tolerance requirements means stopping of the iteration loop, whereas non-fulfilment means that the iteration procedure goes on to the subsequent (k + 1) th iteration.A detailed description of how the above general procedure is managed is given in the following sections.

G
is determined by identifying its distance to the auxiliary surface Σ G in the direction normal to the edge.The applied procedure may be divided into three parts.In the first part, the edge normal vector  In the second part, intersection with the auxiliary surface Σ G is sought in the direction collinear with vector act i k Γ n .As the auxiliary surface Σ G is described discretely by quadrilaterals (Fig. 2), it may be assumed that only one surface element on Σ G is intersected.Being actually related to the edge point act i k P Γ in the corresponding determination of the deviation d i k , let us denote this surface by Σ S i (see Fig. 5) and its nodal points by Σ P j i j ( =1, 2,3, 4) .To enable evaluation of the deviation d i k analytical surface reconstruction of surface Σ S i , Σ F x y z i ( , , ) = 0 , is required, which is carried out in the third part.Apart from the given position vectors Σ P i j j ( =1, 2,3, 4) , the respective surface normal vectors Σ n i j j ( =1, 2,3, 4) at those points are also needed in accordance with Section 2.4.The respective normal vectors' determination follows the procedure described in Section 2.2.
Although the analytical surface reconstruction of surface Σ S i could be done using the global coordinate system ( , , ) x y z , as demonstrated in Section 2.4, it is better to have the surface Σ S i representation defined with respect to a local coordinate system, say ( , , )    x y z , which is built by taking the vector triple act i k After performing all coordinate transformations on the position vectors Σ P i j j ( =1, 2,3, 4) and surface normal vectors Σ n i j j ( =1, 2,3, 4) , as shown in Section 2.2, and obtaining vectors Σ  P i j and Σ  n i j , the geometry of the considered surface element Σ S i can now be functionally approximated in the local coordinate system ( , , )    x y z .By adopting the same functional representation as considered in Section 2.4, we have: As in Eq. ( 12), but expressed in terms of the local coordinate system, the evaluation of the polynomial coefficients  a m m i ( =1, 2,...,12) is based on given topological requirements fulfilment.
Finally, with the interpolation function Σ f x y z i ( , , ) 0    = fully defined, the deviation d i k of the considered edge point act i k P Γ of the actual product geometry G act k can be readily determined.Let us first remember that the coordinate  z axis has been chosen in such a way as to coincide with the edge normal vector act i k Γ n with the origin set at point Γ P i k , i.e. act i k Γ  P = ( , , ) 0 0 0 .Thus, since the deviation is quantified by the distance measured from point act i k P Γ to the auxiliary surface Σ G in the direction normal to the edge, it can be retrieved directly from the functional representation Σ f x y z i ( , , ) 0    = , when taking   x y = = 0 in Eq. ( 13).The established deviation is evidently: and the respective deviation vector d i k is: Fig. 6.Edge geometry deviation determination The deviation of the actual product edge geometry from the reference one, which is established by the deviation vector set d i k { } , determines the appropriateness of the blank shape geometry G bl k used in the forming process simulation in the given iteration.The maximal shape deviation has been chosen as a measure of this appropriateness: The deviation of the actual product edge geometry is computed in each, say, k th iteration.Its determination requires computation of edge normal vectors act i k Γ n at points lying on the edge Γ act k of the actual product geometry G act k and finding a corresponding surface element Σ S i on the auxiliary surface Σ G .The respective surface normal vectors Σ n i j j ( =1, 2,3, 4) at points of the surface element Σ S i , which are also required in computation of the actual product edge geometry deviation, are computed only once, i.e. after the auxiliary surface Σ G determination (see Section 3.1).
A Method for Optimal Blank Shape Determination in Sheet Metal Forming Based on Numerical Simulations

Blank Shape Geometry Correction
In the case where the convergence criterion is not met, i.e. d d max k tol > , a correction of the blank shape geometry G bl k is required.In our approach, this is carried out on two levels separately, first, by carrying out the repositioning of the edge points bl i k P Γ , which is then followed by readjustment of the interior points In order to preserve the quality of the simulation in the subsequent iteration, the aspect ratio of the existing FE mesh, used in the description of the blank shape geometry G bl k , should not be significantly changed by repositioning of the edge points bl i k P Γ .However, this threat can be efficiently neutralised by adequate readjustment of the interior points bl i k P considering the applied edge repositioning, Eq. (17).Numerically, this can be achieved by exposing the blank having shape geometry G to given displacement boundary conditions (Fig. 7a) and solving the corresponding boundary value problem, while assuming isotropic linear elastic behaviour.The blank edge nodes are displaced in accordance with the prescribed correction of the blank shape, Eq. ( 17), while nodes lying on the symmetry planes are subject to symmetric boundary conditions.The result of such computer analysis is the adjusted FE mesh specifying the blank shape geometry G bl k +1 , which will be used in the subsequent iteration.The adjustment is carried out in accordance with the computed displacement field (Fig. 7b).

CASE STUDY -BLANK SHAPE OPTIMISATION
In this section, the developed blank shape optimisation method is validated by considering a particular product with the prescribed edge geometry, manufactured by sheet metal forming (Fig. 8).The product's geometry, which is originally double planar symmetric by design specification, is characterised by two deep cut outs extending downwards to the cup's bottom.Computer simulations supporting the blank shape geometry optimisation in the considered case study have been treated as three-dimensional and carried out using the FEM computer program ABAQUS.The numerical simulation of the complete production cycle of the formed product, which consists of a loading and unloading stage, associated respectively with irreversible plastic deformation and elastic strain relaxation, has been carried out considering particularities in the numerical treatment of the problem.Accordingly, the forming process simulation is carried out using an explicit time integration scheme in ABAQUS Explicit [30], whereas the springback simulation is carried out by ABAQUS Standard [31] using an implicit time integration scheme.For computationally efficient integration of the constitutive equations, the NICE scheme, an explicit integration scheme developed in our laboratory, see Halilovič et al. [3] and Vrh et al. [4], has been implemented into ABAQUS/Explicit via the VUMAT subroutine [32].The optimisation procedure is carried out using a closed loop Fortran program that has been upgraded by an Abaqus' Python script to enable display of the optimisation results in each iteration.
After the optimisation blank shape geometry procedure is completed, which is treated in Section 4.2, both the conceived forming process simulation model and the proposed blank shape geometry optimisation procedure are validated experimentally, which is presented in Section 4.3.

Initial Blank Shape, Geometry and Tool Geometry
The geometry of the considered product (referred as the "reference geometry" in this paper) to be obtained by a sheet metal forming process, where a further sheet metal cutting phase is excluded, is shown in Fig. 8a.An aluminium sheet blank of 1.0 mm thickness is used.The cross-section of the corresponding forming tool assembly is shown in Fig. 8b.The forming tool is defined in such a way that its die geometry is determined by the reference product geometry, with a constant clearance of 1.2 mm ensured between the die and the punch, when the tool is in a closed position.In accordance with the assumptions given in Section 2.5, the tool geometry is considered fixed and not subject to possible variation.The design parameters specifying the initial blank shape geometry, used to start the described numerical optimisation blank shape geometry procedure, are displayed in Fig. 8c.

Determination of Optimal Blank Shape
Determination of the blank shape geometry appropriate for the production of the product, displayed in Fig. 8a, is carried out iteratively following the optimisation procedure, described in Section 3. Thus, a corresponding forming process computer simulation is performed in each iteration, considering the actualised blank shape geometry.
In simulation, the following material and technology process data have been taken into consideration.It is assumed that the product is made from a 1mm thick aluminium sheet exhibiting orthotropic material behaviour, with a Young's modulus of 72500 MPa, a Poisson's ratio of 0. To avoid wrinkling a blank-holder force of 5.2 kN in total is applied.The design prescribed clearance between the die and the punch in a closed position, which is set to 1.2 mm, allows for a potential increase in the sheet thickness during progressing of the punch.For the forming process simulation to be realistic the tribological conditions between the surfaces in contact with one another have to be adequately described.In this investigation the Coulomb friction law is adopted with the friction coefficient between the sheet metal and the tool parts assumed to equal 0.12.The forming process simulations are based on a FE model, which is conceived in such a way as to cope efficiently with all peculiarities encountered when numerically modelling such a complex problem.Since both the tool and the blank shape geometries exhibit planar symmetry with respect to the xz and yz planes and orthotropic material properties are taken into account, only a quarter of the product can be considered in the model (Fig. 10), with symmetry boundary conditions applied and a given blankholder force of magnitude 1.3 kN.In the FE model, the tool surfaces are modelled by rigid FEs with a characteristic size of 1.0 mm, while the sheet metal is modelled by deformable reduced integration quadrilateral shell FEs with a characteristic size of 1.2 mm and with 11 integration points evenly distributed through the shell thickness.
Once the numerical simulation model is completely defined, the steps in the optimisation procedure, given by the flow chart in Fig. 3, can be carried out and the optimal blank shape geometry is achievable iteratively.In Fig. 11, the deviation of the product edge geometry from its reference geometry along the edge contour, as determined by the computer simulation in the considered iteration, is plotted.From these plots it is evident that the developed optimisation procedure is computationally very efficient.Namely, the maximal product edge deviation computed when the initial blank shape is considered in the forming simulation is equal to 3.1 mm.After the first iteration this deviation is reduced to less than 0.3 mm and after two additional iterations to less than 0.001 mm.The shape variation, to which the blank of initial shape geometry is subject when it is optimised using the presented approach, is displayed in Fig. 12.

Experimental Validation of the Blank Shape Optimisation
In order to validate the developed approach to blank shape geometry optimisation in the investigated case study experimentally, the forming tool parts have been manufactured and assembled in accordance with the design scheme (Fig. 8b).In fact, the forming tool is designed in as simple a way as possible.The punch movement is guided by four cylindrical guides, which also have the function of limiting the punch bottom position.The blank-holder is attached to the tool die by eight sets of bolts and springs under tension so that required blank-holder force is ensured.The blanks used in the experiment are cut out of an aluminium sheet by a water jet cutter.
The experiment was carried out on a HI-KON (model HK250S2) forming press with corresponding measurements of the product edge geometry done by the 3D measuring machine Brown & Sharpe, Mistral 7-10-7, its accuracy classified as being ±0.004 mm.The product edge geometry was measured at 6 control points along the edge contour (see Fig. 13 for their position).
To validate the accuracy of the determined blank shape geometry and computed product shape, two experiments were performed, considering different blank shapes.In the first experiment, the product is formed from the blank of initially defined shape geometry (Fig. 8c) in order to confirm the physical and numerical appropriateness of the forming process simulation model definition used in the blank shape optimisation procedure.The difference between experimentally (Fig. 13, Exp.initial) and numerically (Fig. 13, Num.initial) obtained product edge geometry deviation is less than 0.20 mm.The difference between the computer simulation and experiment being in the same range as the accuracy of the cutting machine (±0.1 mm), the adequacy of the simulation model can be confirmed.The second experiment was carried out to confirm the adequacy of the numerically optimised blank shape geometry.As it can be seen from the plots in Fig. 13, the difference between experimentally (Fig. 13, Exp.optimised) and numerically (Fig. 13, Num.optimised) determined product edge geometry is of the same order of magnitude as in the first experiment and is less than 0.22 mm.Therefore, we can confirm that the applied optimisation procedure is appropriate.Photographs of the formed products, produced by the blanks having initial and optimised shape geometry, are displayed in Fig. 14.

CONCLUSION
In this paper, a blank shape optimisation method that enables determination of appropriate blank shape geometry, by which a formed product of prescribed geometry can be produced, has been presented.The developed method is based on an iterative procedure where in each iteration the adequacy of the actual blank shape geometry is checked by identifying the edge geometry deviation of the resulting product geometry, which is obtained by means of the

1 e and v 2 eFig. 1 .
Fig. 1.Surface normal vector determination: a) FE's nodal normal vector n i e determination, and b) auxiliary vector ′′ n i determination

i i ≡ 1 (
Fig. 1b).The base vectors e e e defining the respective local coordinate system ( , , )    x y z at point P i in accordance with the above stated properties, are determined by expressions: e e n r e e e e n , r x , r y and r z are the components of vector r = ( , , ) r r rx y z and quantity v( ) θ is defined by v( ) = 1 ( ) θ θ − cos , the angle θ denoting the angle between the base vectors e z and e  z .To enable analytical surface S i representation f x y z i ( , , ) = 0    with respect to the local coordinate system ( , , )

Fig. 2 .
Fig. 2. Edge normal vectorΓ n i determination Σ G The appropriateness of the actual blank shape geometry G bl k can be established by measuring the deviation of the actual product edge geometry Γ act k from the reference edge geometry Γ ref , which can be done in many ways.Here, the deviation will not be measured directly with respect to Γ ref , but indirectly by making use of an auxiliary virtual surface Σ G emanating from the reference product edge Γ ref .This surface can be generated considering surface topological properties of the reference product geometry G ref on its edge Γ ref .Lines l i , aligned along the surface normal vectors ref i n at points ref i P Γ of the discretised edge, are used, accordingly, as the building basis for the auxiliary surface construction (Fig. 4).To enable numerical control of the actual product edge adequacy this surface is further discretised by quadrilateral sub-domains, yielding point topology Σ G .Having those properties, i.e. emanating from points ref i P Γ of the discretised reference product edge Γ ref ( Γ ref ⊂ Σ G ) and being at those points perpendicular to the reference product geometry G ref , the surface Σ G may be considered in the iteration process as the target surface, which is to be attained by the edge points act i k P Γ of the actual product G act k , as closely as possible.For a determination of the auxiliary surface Σ G , which is computed before starting the optimisation procedure, computation of the surface normal vectors at points lying on the edge Γ ref of the reference product geometry G ref is required.

Fig. 4 .
Fig. 4. Auxiliary surface Σ G determination accordance with the procedure, described in Section 2.3.With the surface normal vector act i k n and tangential vector act i k Γ t computed by Eqs.(8) and (9), respectively, considering surface and edge topology properties of the actual product geometry G act k , the vector act i k Γ n is determined by Eq. (10).

Fig. 5 .
Fig. 5. Finding surface element Σ S i on Σ G intersected by a line collinear with vector act i k Γ n . In this way, occurrence of round-off error due to possible computing with large numbers is avoided and accuracy of the computation ensured.
bl i k P appertaining to the FE mesh.Correction of the edge blank shape geometry Γ G bl k is carried out based on the established deviation vector set d i k { } , Eq. (15), with repositioning of points bl i k P Γ applied in the direction of the edge normal vectors bl i k Γ n , Eq. (10), considering deviation magnitudes d i k .The positions of the points lying on the edge Γ bl k +1 of the newly defined blank shape are thus computed by the equation:

Fig. 7 .
Fig. 7. Blank shape geometry correction G G bl k bl k → +1 ; a) applied displacement boundary conditions, and b) adjustment of the FE's mesh

Fig. 8 .
3, an initial yield limit of 105 MPa, and coefficients describing orthotropic material properties in relation to the potential function Hill 48 being: R xx = 1 0 ., R yy = 0 95 .and R xy = 0 97 . .The yield curve specifying the hardening plastic behaviour of the given material is plotted in Fig. 9. Since no significant reversed plastic bending occurs during the considered forming process, possible kinematic hardening is neglected and isotropic hardening is assumed.Product and forming process geometry design specification; a) product reference geometry, b) forming tool assembly, and c) initial blank shape geometry