Optimization of the countersinking parameters based on the response surface method

The countersinking process is affected by many factors including the tools and the workpiece parameters. Some forming phenomena such as the knife-edge affect the quality of the countersunk hole. Up to now, many kinds of research rely mainly on experiments which lead to poor quality and difficult control of this process. In this paper, a proposed numerical optimization of the countersinking process is developed to obtain a normalized countersunk hole. This optimization approach is based on the response surface method (RSM), design of experiments (DOE) and sequential quadratic programming (SQP). The finite element model is performed with an elasto-plastic behaviour for simulating the process. A configuration with an imposed displacement applied to the blank holder is adopted in this study. The comparison between the results of the numerical model and the experiments showed a good agreement.


Introduction
In the industries, such as aircraft manufactories, the countersunk hole is widely used to joint plates with rivets or screws.
The countersinking process using forming techniques improves productivity and decreases the cost and gives good mechanical properties with no waste of material. Therefore, the quality of the countersunk hole plays an interesting role to ensure the rigidity and the accuracy of the joined plates.
Many recent pieces of research have been focused on crack behaviour and fatigue life. [1] investigated the cracking behaviour of the fastened plates using rivet and countersunk holes. Low and high loading cases were discussed. They found that the beginning of the crack and its location in the countersunk hole depend on the degree of the load transfer through the fastener in the join specimens. [2] investigated the effect of the detail design on the fatigue life of the countersunk hole in the joining sheet metal plates. They found that the fastening using countersunk hole increases the fatigue life of the sheet assembly more than the use of the traditional drilled hole.
In other works, [3] studied the effect of the shape of the countersunk hole on the fatigue life of the fastened sheet metal plates. A comparative study was performed between different geometries of the countersunk hole. The first one has a rounded corner shape and the second was with a simple shape. It was found that the fatigue life was enhanced by using the rounded corner. Some other works are interested in studying the stress concentration factors in a countersunk hole in plates [4]. In their work, they carried out a parametric analysis to study the main effects of some parameters on the stress concentration. As consequence, they determined a new equation to model the factor of this concentration of the stress around the countersunk hole. In addition, [5] studied the effect of the surface of the countersunk hole after the manufacturing process. They found that metal protuberances are present in the manufactured mini-plate which are used in bone operation surgery and they cause inflammation.
Hassen et al. [6] studied the forming kinematic during the countersinking process and they defined many new phenomena such as the expansions and the gaps. The effect of the pre-hole fillet is investigated too. They found that the use of an adequate fillet in the numerical model leads to obtain acceptable results comparing to the experimental ones.
Hassen et al. [7] studied the effects of the geometrical parameters of the pre-hole on the final shape of the obtained countersunk hole. They used a statistical technic based on the definitive screening design and the diagram of the main effects to perform a sensitivity analysis. In this work, they carried out a parametric simulation analysis using three parameters: the radius of the pre-hole r i , the fillet radius r t of it and the thickness of the workpiece e. They found that the r i has the dominant effect on the final shape of the countersunk hole as well as on the punch load. The parameter r t has a notable significant effect on the punch load. Concerning the thickness e, it has a significant effect just on the height h c of the countersunk hole. They plotted the curves of many forming kinematic with variable values of the studied parameters to better understand their evolutions during the countersinking process and to confirm the statistical results. In their works, they carried out experiments to validate their numerical results.
Up to now, there are few works about optimization in the countersinking process. In [8], a graphical optimization method is presented to obtain standard countersunk holes. They gathered the results of a parametric study in the same graph. By using a complex projection, the punch travel and the initial hole diameter can be extracted from the graph to obtain the wished normalized countersink. This graphical method could not arrange all the forming kinematic phenomena and leads to losing the accuracy of the results because of the complexity of the layout of the curves in the same graph and the projection method.
Numerical optimization is a powerful technique, which is largely used in many fields of researches. In the forming process, the cost functions and the constraints of the optimization problem are rarely given with analytical form which is always hard to obtain. The response surface methodology (RSM) is used to overcome these problems by developing approximate functions based on results calculated with parametric numerical simulations [9]. The numerical design of experiments (DOE) is used to organize a set of points of the design space for which the numerical simulations are performed. When using the RSM method, the best choice of the DOE and the function of approximation lead to having accurate results [10]. The validation of the calculated response is usually confirmed by statistical coefficients [11] [12].
In this paper, the effects of some factors on the forming kinematics and the loads are numerically investigated. The adopted factors are the punch travel δ p and the initial inner radius r i of the workpiece. This paper aims to search for the optimal values of these factors to obtain a normalized countersunk hole. The developed optimization approach is based on the RSM methodology coupled with the SQP algorithm to resolve the complex cost functions. The used approximate functions are obtained by the RSM technic based on numerical simulations according to definitive screening design. The cost function is defined as the quadratic error between the calculated values and the average normalized parameters of the countersink geometry. The constraint function is defined to prevent the edge-knife phenomena.
The adopted approach showed good accuracy and reliability to obtain optimal solutions. The numerical results are validated by experiments.

Finite element modelling and parameters
The studied model was carried out in Abaqus/Standard®. To establish this model, it is necessary to define the material behaviour, the parameters of the parts and the boundary conditions.

Material behaviour
The material used in this study is the steel S235. The mechanical properties of the metal are listed in Table 1. Since the countersinking process could be modelled as a compression state, the uniaxial compression test was used to describe the stress-strain relationship. This test is carried out by the use of the technique described in recent works [6,13,14]. The compression test was performed using specimens stack with 5 discs. All the discs were cut from the adopted sheet metal with a thickness of 2.9mm and a diameter equal to 14.5mm. Therefore, the obtained specimens were similar to a cylinder with a diameter and a length equal to 14.5 mm. The obtained stress-strain curve is shown in Fig. 1. The used material was assumed to have an elasto-plastic behaviour and to obey the von Mises yield criteria with isotropic hardening rule. For the adopted FE model, a multilinear uniaxial stress-strain was used in Abaqus®. (Fig. 2)

Countersinking parameters
The used workpiece is a washer cut from steel (S235) sheet metal (Fig. 3). The thickness e=2.7mm is considered for the washer. Its external radius r e is fixed equal to 15mm and the internal radius r i is assumed to be calculated using the adopted optimization method. The fillet radius of the pre-hole is set to be r t =0.1mm. The entire parameters of the process are listed in Tables 2 and 3.
During the countersinking process, this initial workpiece has followed many kinds of forming kinematics. Therefore, new parameters were defined to describe and quantify the observed phenomenon in the final shape of the workpiece. The entire parameters are presented in Fig. 3 and listed in Table 3.
In Fig. 3, the punch was characterized by its radius r p and its conical angel Φ p . The blank holder has an internal radius r b =5mm with a fillet radius equal to r s =0.1mm. All the dimensions of the tools are detailed in Table 2.

Finite element model and boundary conditions
In the present work, the finite element model is performed as an axisymmetric model as shown in Fig. 2. The workpiece is modelled as a deformable body. The punch, the blank holder and the die are modelled as rigid bodies. The die is clamped and the punch is constrained to have only an axial move. Before the punch goes through the inner hole of the workpiece, an imposed displacement was applied on the blank holder to prevent the upward movement of the workpiece during the countersinking process. To enhance the countersinking process, a non-null displacement (δ b =1μm) is applied on the blank holder instead of the null value (δ b =0μm) which was used in the maintained configuration of the previous work [6,7].

Meshing
Since the FE model is considered an axisymmetric assembly, the used mesh is assumed to be axisymmetric. To reduce the time of simulations, the mesh is chosen as CAX3H with reduced integration. Due to the large strain level in countersinking process and to avoid the distortion phenomenon, the size of the mesh is resolved by using the adaptive process in Abaqus/Standard®. An adaptive remeshing rule, based on the error indicator, is used to have an optimal mesh after 3 iterations. A typical mesh is presented in Fig. 2.
The contacts between the tools and the workpiece are modelled with Coulomb's friction law. To describe reasonable contact conditions, the friction coefficient is fixed equal to μ=0.12 [6,7].

Validation of the finite element model
The described finite element model is validated by comparing the results of the experiments and the simulations.
The experiments are carried out using a workpiece with a thickness e=2.7mm and an external radius r e =15mm. The prehole is drilled to have an inner radius equal to r i =2.8mm. The performed experiments are used to validate the new FE model in which the imposed displacement on the blank holder was non-null.
The experimental tools are manufactured by respecting the FE model parts as shown in Fig. 4. The experiment process is detailed as the following: -The experimental tests were carried out with a universal tension-compression testing machine of maximum load 200 kN. -All the experiments were performed with a controlled displacement and a fixed velocity equal to 5mm/min. -No lubricating was used.
-Firstly, the workpiece was carefully placed in such a way to be coaxial with the punch, the blank holder and the die. -An imposed force equal to 50 KN was applied on the blank holder by using the screw-nut system (screw is the die and the nut is the blank holder) of the experimental tools. -The tools experiments were placed coaxial to the axil of the punch of the compression testing machine. After each countersinking operation, the obtained results were saved in a text file.

Formulation of the optimization problem
As known, the goal of the optimization process is usually to determine the optimum which usually can be the minimum or the maximum of a cost function depending on the design variables [9] (Table 4).   In this works, the objective is to determine the optimal values of the initial radius of the workpiece r i and the punch travel δ p to obtain the normalized countersunk hole. This objective could be modelled by minimizing the quadratic error between both calculated and normalized (average value) values of the countersunk hole parameters as given in Table 5. Besides, the knife-edge is a phenomenon observed when the countersink height exceeds the skin thickness of the plate [15]. Generally, this phenomenon leads to a high stress concentration factor in the countersunk hole. In practice, the height of the knife-edge must be superior to the 1/5 of the plate thickness to avoid the excessive stress concentration and the crack phenomenon [16].
Therefore, the proposed problem can be stated mathematically as follows: & e f is the approximate function of the quadratic error. & δ p and r i are the adopted design variables (N=2 variables). The intervals of these two variables are defined in Table 4. & e r c e r fi and e h c represent respectively the approximate surface response of the superior radius r c , the final inner radius r fi and the height h c of the countersinking in terms of the design variables (δ and r i ).
Concerning the approximate surface responses, they are obtained by the use of the RSM methodology based on the numerical design of experiments DOE and the second-order regression.

Response surface methodology (RSM)
Response surface methodology is a statistical technique that is useful in applications where the responses of interest are affected by several variables. The aim objective of this method is to determine the approximate functions of these responses instead of complex analytical functions.
In this present study, the RSM is used to obtain the approximate responses y i in terms of the design variables (δ p and r i ) as given in Table 4. Generally, the order of the used polynomial approximation is primordial to obtain the best approximate surfaces which fit the simulation data of the studied responses y i . In [17], the second-order polynomial approximation is adopted to predict these responses y i . This model is widely used because of its flexibility to predict the functions of the surface of responses in the small region of the variables space. This mathematical model is sufficient for modelling the curvature region in the surfaces of responses. It is given as follows: where x i is the design variables, β i is the regression coefficients and ε is the approximation error.
In the case of two variables, Eq. (2) is presented as follows: By changing x 3 ¼ x 2 1 , x 4 ¼ x 2 2 and x 5 = x 1 x 2 , the recent Eq. (2) can be converted to a linear regression model as follows: Now, the regression coefficients can be estimated by using the least-squares method as explained in [11,12,17].
The statistical test is usually used to validate the fitness of the calculated response surfaces model. The analysis of variance with ANOVA is one of the powerful tools for this kind of tests. The well-known statistical parameters for evaluating the fitness of the model are R 2 , and R 2 pred . These parameters are explained and detailed in [11,12]. Generally, the lager value of these three statistical parameters, the better fit of the model. Therefore, the approximate model is more accurate when these statistical parameters are closer to 1. In addition, the parameter R 2 increases when the number of terms in the model increases. On the contrary, R 2 adj decreases when the number of unnecessary terms in the model increases. Also, the parameter R 2 pred is used to evaluate the predictability of the model. In the present work, these three parameters are adopted to validate the fitness of all the studied approximate responses y i .

Design of experiments (DOE)
The previous approximations responses y i are based on numerical simulations which are defined in the space of the adopted design variables (δ p and r i ) as presented in Table 4. The DOE is built with a set of points to get the good fitness of the approximate surface responses. In the literature, it could be found many types of DOE such as full factorial, central composite, Box-Behnken and so on. Generally, the Latin hypercube design (LHD) of experiments is widely used for fitting a second-order response surface model with high accuracy [18]. In this work, the used LHD is built with 9 experiments as presented in Table 6. The criteria of maximizing the inter-site distances are used to obtain the optimal LHD as shown in Fig. 5. 6 Results and discussions

Load and forming kinematics
The described FE model is validated by comparison to the experimental results. In Fig. 6, a comparison between the experiment and the numerical final shape of the workpiece is presented. It is noticed clearly that the two results are in good agreement. In Fig. 7, the experimental and the numerical punch force evolution could be considered in good agreement too despite the difference between them in the last stage of the process. The evolution of this load is approximately the same in the beginning and the second stage of the process. In the final stage, a deviation between these two evolutions was observed and a difference about 14% is reached in the end of the process. This difference could be explained by the experimental contact friction which is affected by the contact pressure, the local stress state, the roughness of the contact surface and the sliding velocity. All these factors increase the friction and its effect in the last stage of the experimental process. Contrary, in the adopted FE model, the friction is assumed to be constant, to simplify the model without significant effects on the final countersunk hole shape and the punch load.
Due to the complex flow of the metal between the tools (punch, blank holder and die), much of forming kinematics are observed clearly in the numerical simulations. A descriptive final shape of the deformed workpiece is shown in Fig. 3. To investigate these phenomena, a numerical simulation is performed using the same workpiece used in the previous experiments.
The obtained numerical shape of the countersunk hole is shown in Fig. 8. As can be seen from this figure, the main displacement is located under the punch with high concentrated stress. This local deformation leads to obtain the shape of the countersunk hole which is defined by three important parameters: the final radius r fi , the superior radius of countersink r c and the height h c .
On the other hand, it is noticed the presence of axial expansion just close the countersunk hole, in the upper surface of the deformed workpiece as shown in Fig. 3. This axial expansion is quantified by g h . Radial expansion is observed in the outer radial surface of the workpiece. It is defined as the difference between the final and the initial external radius (r fe -r e ). These two expansions are considered geometric defaults in the countersinking process because they affect negatively the quality of the final shape.
As consequence, it is worth adopting the developed FE model to predict the final shape of the workpiece and the force of the punch during the countersinking process. This potential  5 Design of experiment type LHD for δ p and r i advantage is the key to carry out the RSM method and perform numerical optimization. The shape of the countersunk hole and the punch load will be used for the rest of this study.

Responses surfaces
In the first step, the simulations are executed using Abaqus/ Standard®. The scripting language of Python® is used to extract the responses data y i after finishing all executions. These gathered data are extracted from the last frame of the countersinking process. The obtained results of each response y i are presented in Table 6.
Furthermore, the approximate surfaces responses for each response y i are presented in Fig. 9 and the coefficients of regression are listed in Table 7. The statistical parameters show good quality for all the approximations. Concerning the response h c , the values of R 2 , R 2 adj andR 2 pred are slightly low compared to the other responses, but its approximate function is good enough for the rest of the study. We can conclude that the adopted DOE and the quadratic regression model leaded to satisfying accurate approximate functions for all the studied responses. These approximations will be used in the optimization process in the next section.

Optimization
Firstly, it is worth noting that the response surface plots in Fig. 9 could be used to analyse the effects of the adopted parameters on the obtained countersunk hole and the developed punch force. In this figure, the evolution for each response surface has a slight quadratic zone. Only in the case of the radial expansion g h , high non-linearity is observed and this could explain the complexity of the forming kinematic in this zone.
Concerning the punch force, it is mainly affected by the punch travel δ p more than the radius of the pre-hole r i . This force increases by increasing these two parameters Fig. 9a. The increase of the punch travel increases the surface contact between the punch and the deformed pre-hole of the washer. Since the punch force depends on the surface contact and the applied pressure, this force is increased quasilinearly as shown in Fig. 7. It is worth noticing that    maximum punch force is obtained when using a high value of r i and δ p and vice versa.
Since the shape of the countersunk hole is defined by the responses r fi , h c and r c , the effect of r i and δ p on these responses must be investigated. The punch travel δ p has the dominant effects on h c and r c as shown in Fig. 9c and d. The increasing of this parameter increases linearly these two responses [6,7]. The parameter r i has just only a significant effect on r c . On the other hand, the parameter r i has a dominant effect on the retraction r fi . It is obvious that with the increase of r i , the final value of r fi increases too (Fig. 9b). Concerning the radial expansion r fe -r e , it is mainly affected only by δ p . This phenomenon was observed on the final stage of the process and its value increases notably in the end of the punching travel [6,7].
As consequence, this sensitivity analysis showed the significant effects of these two adopted parameters (r i and δ p ) on the countersinking process and exactly on the final shape of the countersunk hole. This analysis confirms the use of these parameters in the next task of the optimization process as well as in the recent work of Jallouli et al. [8].
After this sensitivity analysis, a proposed optimization methodology is applied for the normalized large series of the countersunk screw M4 as described in norm NF EN 20273 [8 Jallouli]. This methodology is based on the following steps: -Since the approximate functions of the responses y i are determined by the RSM method, the cost function could be defined as described in Equation (1).
-The SQP algorithm has good capability to converge so fast into the global optimum solution. For this reason, it is used to resolve the studied cost function. -The optimum solution is obtained after a few iterations and the constrained function is respected without violation. The obtained optimal design variables are δ p =1.75mm and r i =2.64mm.
These optimal parameters are integrated into the FE model. The obtained countersunk hole is showed in Fig. 10. The results of the simulation show good agreement with the approximate results as shown in Table 8. These results confirm that the geometrical parameters of the optimal countersunk hole are in the normalized interval. Therefore, the developed methodology could be used explicitly to predict the optimal parameters δ p and r i to obtain the normalized countersunk hole.

Conclusion
This study was an attempt to determine the optimal radius of the pre-hole r i of the workpiece and the punch travel δ p to obtain a normalized countersunk hole. An optimization methodology was developed based on the RSM method, the Latin hypercube DOE and the SQP algorithm. Based on the DOE, parametric simulations are performed to determine approximately the studied responses y i using the RSM method. For the studied optimization problem, the cost function is defined to minimize the quadratic error between the desired and the calculated normalized countersunk geometry in terms of r i and δ p . By using the SQP algorithm, the optimization problem is resolved under the constraint of preventing the knife-edge phenomenon.
Some of the following interesting results are obtained: -The adopted FE model with an imposed displacement on the blank holder is worth to be used to predict the final shape and the developed forces for the countersinking process. This FE model is validated by comparison with the experimental results. -The use of the RSM method and the Latin hypercube DOE seems to be a good alternative to obtain an accurate approximation for the studied responses y i . -The plot of the surface responses could explain the impact of the design variables on each response y i . It was found that the punch travel δ p is the dominant factor that affects all the responses. On other hand, the initial radius of the pre-hole r i seems to have a dominant effect on the retraction of the pre-hole r fi and a significant impact on the superior radius of the countersunk hole r c . The effect of r i becomes less important on the punch force and the expansions compared to the effect of the punch travel δ p . -The developed optimization approach provides the designer with accurate results for the normalized countersunk hole.
The future objectives consist of improving the developed optimization method. It still needs the use of more parameters to obtain a normalized countersunk hole for all the series (large, medium and small) under other technology constraints and damage limitations. Besides, the use of the anisotropy of the material will be a good step to enhance the FE model. The improvement of the adequate friction law will be another key to enhance the FE model results.
Author contribution This manuscript is truthful original work.

Declarations
Consent to participate Informed consent was obtained from all individual participants included in this study.
Consent for publication All individual participants declare that they accept to submit this paper for publication.

Conflicts of interest
The authors declare no competing interests.