**Journal of Materials Science and Chemical Engineering**

Vol.04 No.04(2016), Article ID:66113,11 pages

10.4236/msce.2016.44002

Hardness Profile Prediction for a 4340 Steel Spline Shaft Heat Treated by Laser Using a 3D Modeling and Experimental Validation

Mahdi Hadhri, Abderazzak El Ouafi^{*}, Noureddine Barka

Computer Science and Engineering Department, University of Quebec, Rimouski, Canada

Copyright © 2016 by authors and Scientific Research Publishing Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY).

Received 6 April 2016; accepted 25 April 2016; published 28 April 2016

ABSTRACT

Laser surface transformation hardening becomes one of the most effective processes used to improve wear and fatigue resistance of mechanical parts. In this process, the material physicochemical properties and the heating system parameters have significant effects on the characteristics of the hardened surface. To appropriately exploit the benefits presented by the laser surface hardening, it is necessary to develop a comprehensive strategy to control the process variables in order to produce desired hardened surface attributes without being forced to use the traditional and fastidious trial and error procedures. The paper presents a study of hardness profile predictive modeling and experimental validation for spline shafts using a 3D model. The proposed approach is based on thermal and metallurgical simulations, experimental investigations and statistical analysis to build the prediction model. The simulation of the hardening process is carried out using 3D finite element model on commercial software. The model is used to estimate the temperature distribution and the hardness profile attributes for various hardening parameters, such as laser power, shaft rotation speed and scanning speed. The experimental calibration and validation of the model are performed on a 3 kW Nd:Yag laser system using a structured experimental design and confirmed statistical analysis tools. The results reveal that the model can provide not only a consistent and accurate prediction of temperature distribution and hardness profile characteristics under variable hardening parameters and conditions but also a comprehensive and quantitative analysis of process parameters effects. The modelling results show a great concordance between predicted and measured values for the dimensions of hardened zones.

**Keywords:**

Heat Treatment, Laser Surface Transformation Hardening, Finite Element Method, Hardness Profile Prediction, AISI 4340, Nd:Yag Laser System, ANOVA

1. Introduction

Surface transformation hardening processes are designed to improve wear and fatigue resistance by hardening the superficial critical areas using brief and localized heat gains. Among these processes, laser surface transformation hardening process is well-known by his capacity in terms of power flux density and is recognized by his fast, local and accurate thermal cycles, while limiting the risks of undesirable distortion and deformation effects. Application of the laser beam rapidly raises the surface temperature (more than to 1000 K/s), resulting in a thin layer that is converted into austenite. Subsequent removal of this energy results in self-quenching caused by the conduction of heat into the relatively cool bulk of the material. This produces a rapidly cooled surface layer and causes a transformation of the austenite into martensite [1] - [5] .

Laser surface transformation hardening offers several advantages: localized treated areas, a relatively small heat affected zone, limited metallurgical changes, reduced residual stresses, very fast thermal cycling and autogenous quenching, and appropriate process for automation and complex production lines when using robots. Despite all its advantages, applications of laser surface treatments represent a very small percentage of industrial plants. Laser surface hardening is still in its infancy, with only a few years of development and it is virtually only developed in the aerospace and automotive industries. This work is a continuing effort to develop power laser applications for surface treatments.

The modeling of laser heat treatment went through several stages. In the first, researchers were interested in statistical modeling to understand the influence of certain parameters and develop empirical formulas [6] - [8] . Then, other researchers became interested in analytical modeling based on the general equation of heat conduction proposed by Fourier [9] - [11] . During the 90s, advancements in IT brought more powerful computational tools to researchers in all fields. Several numerical modeling platforms made their appearances and greatly accelerated technical developments. These improvements allowed more complex problems to be modeled and solved. According to the literature, there were three methods to model a mobile heat source. The first method was based on the Rosenthal equation of a mobile heat source [12] [13] . The second was to move the source implicitly based on the transport term in the heat equation [14] . The last one was the method of Area Sector Approach. The geometry of the heat treated sample was the main factor in selecting the method of modeling [15] [16] . The second method was efficient only if the geometry was not intricate, and the third was applicable when the movement was simple. The first method was effectual with complex geometries and complex movement, but the challenge was finding the ideal trajectories of the heat source.

To treat revolutionary geometries, the sample must be rotated, and the laser beam must be in a translational movement as illustrated in Figure 1. In accordance with the literature, Rahul Patwa and Shin achieved a 3D finite-element model [17] . The model combines a transient digital three-dimensional solution (based on the modeling of Rozzi et al. [18] [19] ) for a rotary cylinder undergoing laser heating by beam translation with a kinetic model. In order to verify the results of the simulation, an experiment is performed. Both researchers reach a depth of 0.54 mm with a hardness of 63 HRC on an AISI 5150 steel sample with a laser (diode) power of 500 W and a rotational speed of 6 RPM [17] . Skvarenina et al. were capable of predicting and experimentally validating a 2.5 mm hardening depth with a uniform hardness of 57 HRC on an AISI 1536 steel cylinder 60 mm in diameter, using a scanning speed of 2.9 mm/s, a diode laser power of 1220 W and a rotation speed of 1 RPM [20] . Another thermal transient 3D model is developed by Leonardo Orazi et al. [21] . The model is based on the geometry of the ring spot and was validated by experimental tests. The advantage of the Leonardo model over other models is that it achieves very high speeds. For a rotational speed of 1140 RPM, a power of 1 kW, a scanning speed of 30 mm/min, and a test piece of AISI 1040 steel 30 mm in diameter, he found a hardness of 690 HV. In general, a second laser pass generates a tempering of the material that is characterized by a drop in micro hardness. In the same context, low processing speeds create a superposition of treatment which gives a non- homogeneous micro hardness.

The literature review reveals the small number of researches dealing with revolutionarily complex sample processing. The majority of the researchers focused mainly on the study of this phenomenon on gears. Benedict and Eskildsen tested an approach to treat small gears, which proved very promising [22] . In fact, this method consisted of the laser beam scanning from one gear tooth tip to another through variation of the angle of incidence, power, and interaction time (forward speed). In 2003, Zhang et al. used this approach to treat sprockets with 98 mm outside diameters and 23 teeth [23] . They were installed on a mounting allowing the wheel to be moved laterally and be brought into rotation. The results were very conclusive: the cured depth of the flanks was

Figure 1. Basic configuration of the laser hardening of a cylindrical part.

of a relatively uniform thickness and was not a stressed fusion surface. However, no process modeling was presented. Pretorius and Vollertsen have modeled a form of laser heat treatment on a toothed wheel, where the heat treatment was applied only in the tooth root [24] . The 3D digital model was developed using the SYSWELD package and consisted of modeling the thermal flow, metallurgical transformations, and geometrical. Clearly, the modeling of heat treatment of the revolutionary complex geometry is limited. Among such geometries are spline shafts. Until now, there is no work, neither experimental nor modeling, which presents the laser heat treatment of the spline shaft. Due to its ability to transmit large torques and ease machining, spline shafts have become the essential tools for power transmission. In this work, a method for the prediction of hardened depth using laser heat treatment of a spline shaft with a high speed of revolution is presented.

The main objective of this work is to develop an integrated approach for hardness profile predictive modeling and experimental validation for spline shafts using a 3D model. The numerical simulation of the hardening process is carried out by 3D finite element model using Comsol Multiphysics software. The model is used to estimate the temperature distribution and the hardness profile attributes for various hardening parameters and material properties. Applied on AISI 4340 steel spline shaft, the experimental calibration and validation of the model is performed on a commercial 3 kW Nd:Yag laser system using a structured experimental design and confirmed statistical analysis tools. The results reveal that the model can provide not only a consistent and accurate prediction of temperature distribution and hardness profile characteristics under variable hardening parameters and conditions but also a comprehensive analysis of process parameters effects. The results show great concordance between predicted and measured values for the dimensions of hardened zones.

2. Finite Element Modeling

In this study, the laser was modeled as a source of circular Gaussian heat. The laser moves along the spline shaft. The latter is mounted in a test stand that allows it to be turned at rather high rotational speeds.

In general, non-linear mathematical models of heat transfer by conduction in a homogeneous and isotropic medium take the following form:

(1)

The volume density of the laser Q (x, y, z, t) applied to the material, is given by:

(2)

Here Q_{0} is the power of laser, A_{c} is the coefficient of absorption, R_{c} is the coefficient of reflection, W is the radius of the laser beam, and f (x, y, z, t) is the function that describes the shape and the path of the beam, given by:

(3)

The challenge in modeling the laser heat treatment of splines is that the geometry is complex and the heat source must follow the teeth and the flanks of groove. For this, the solution was to find a mobile frame of vector space at each time t, of which a, b and c are given by:

(4)

where x_{1} is the length travelled by the laser beam, t_{f} is the duration of treatment, ω is the rotational speed in rad/s, y_{0} and z_{0} are the instantaneous positions along y and z, respectively, and are given by Equation (5):

(5)

The beam should treat tooth and flank, which requires clarified beam position either on the tooth or the flank.

(6)

Here R_{d} is the outer diameter of spline shaft, e is the depth of the groove and H (t) is the position function of the laser beam:

(7)

where A is a function dependent on the number of teeth.

3. Metallurgical Modeling

The metallurgical transformation process for the heat treatment of steel occurs over three major steps: the pearlite transformation to austenite (pearlite dissolution), the homogenization of the carbon in austenite, and the austenite transformation to martensitic [25] . By heating the material up to the temperature of eutectoid Ac1, colonies of pearlite in the microstructure are transformed into austenite. The distance between the pearlite plates, which allow colonies pearlite to be completely transformed into austenite, is given by the following formula:

(8)

where D_{0} is the diffusion constant, Q_{a} is the activation energy, R the gas constant, T_{p} is the temperature of spades, and the two constants c_{d} and are given by Equations (9) and (10):

(9)

(10)

Here K the thermal conductivity, V is the scanning speed, T_{0} is the initial temperature and R is the reflection coefficient. The homogenization mechanism is simple: around a ferrite grain and a cementite grain, an austenite germ can be created. This germ is formed by eutectoid transformation with a chemical composition of 0.8% C. As the temperature rises it undergoes a systematic change in its composition. Rapid cooling of the austenite, which is formed only within a thin layer during laser hardening due to the self-sealing of the material when the laser beam is moved away, makes it difficult for carbon to diffuse outside its lattice [26] [27] . When the carbon is trapped in the network and cooled, the face-centered cubic crystal structure of austenite is transformed into a hybrid quadratic structure, called martensite [5] . The martensitic volume fraction, f, which is formed on a period T, is given by Equation (11):

(11)

where f_{i} = C/0.8 as the initial volume fraction of pearlite and f_{m} is the volume fraction of martensite, given by the following relationship:

(12)

The hardness of the material is calculated as follows:

(13)

H_{m} and H_{f+p} are calculated following Maynier’s equations and taking into account the initial chemical composition [10] [11] .

4. Simulation Results

This study investigates the machine sensitivity parameters of laser heat treatment of an AISI 4340 steel spline shaft. The AISI 4340 steel is very common in the aerospace and automotive industries in the manufacture of propeller shafts, connecting rods, gear shafts and other parts, and automobiles due to its high tensile strength. The AISI 4340 chemical composition is given in Table 1. The sample used in the simulations and validation is 15 mm in diameter, 2.5 mm in thickness and is inclined at an angle of 20˚C. Primary results of the simulation of heat treatment showed that there is no heat affected area and the temperature does not exceed 350˚C, even with a speed of 2 mm/s, a power of 2500 W, and rotation speed of 1500 RPM. The results are experimentally validated. To achieve a heat affected zone by the laser beam, it must pass through a preliminary heating. Multiple scanning was carried out to increase the original sample temperature from 20˚C to 500˚C. Table 2 shows the material properties used for the simulation.

It is clear in Figure 2 that the laser is in the process of turning around the spline shaft, treating both the tooth and the flank. The temperature rises to 837˚C with a rotational speed of 1000 RPM, a scanning speed of 5 mm/s, and a power of 2200 W. Note that the austenitizing temperature is 790˚C.

Figure 3 and Figure 4 show that the temperature progressively increases approaching the measurement point,

Table 1. Material chemical composition of 4340 steel.

Figure 2. Distribution of temperature.

Figure 3. Distribution of heat flux.

Figure 4. The evolution of temperature versus time (x = 5, y = 0, z = 15).

then gradually decreases away from the measurement point. It is clear that the temperature increases and decreases when the laser beam turns around the spline shaft.

5. Experimental Validation

A Yag laser with a maximum power of 3 KW was used to validate the finite element model. The laser head is mounted on a Fanuc robot with six degrees of freedom (see Figure 5 left). The laser beam diameter is evaluated at 1.08 mm when focused. The specimens are treated by a hardening and tempering process to ensure a core hardness of 35 HRC. The latter is the untreated area. The prediction algorithm is intended to look at the hardened area. Two validation tests shown in Table 3, are performed to validate the hardened depth at the tooth (D_{t}) and the flank (D_{f}).

The curve describing the hardness profile is divided into three zones. The first is the hardened zone, consisting of 100% martensite. The second is the transition zone, consisting of ferrite, perlite and martensite. The third

Table 2. Material properties of 4340 steel.

Table 3. Numerical model validation tests.

Robotic laser cell Sample of hardness profile-2200 W, 4 mm/s and 2000 RPM

Figure 5. Experimental setup and example of hardness profile.

zone is the untreated area. The transition zone is not considered in the modeling. The modeling of the transition zone with precision depends on two main parameters: the cooling rate and the initial hardness of the material [26] [27] . In the Figures 5-7, we can see that the error of modeling the martensitic area is below 15%. The tests show that the flank hardened depth is always less than or equal to the depth of the hardened tooth.

6. Statistical Analysis

The objective of this part is to identify the influence of various system parameters on the hardened depth. This is carried out with the help of experiment designs, consisting in producing a series of N experiments and determining the value of the response function for these N configurations. So, in this case, the selected solution is the profile of the hardened depth shown by the tow characteristic measurements {D_{t}, D_{f}}. The experiments are carried

Figure 6. Sample of hardness profile-2200 W, 4 mm/s and 2000 RPM.

Figure 7. Sample of hardness profile-2200 W, 4 mm/s and 2000 RPM.

out using a 3D model with COMSOL software (see Figure 2). The factors to be examined in this study are power, rotary speed and scanning speed. Taguchi proposed that in order to optimize a process or a product, experiments should be carried out in a three-step approach, i.e. system design, parameter design, and tolerance design [28] . The Taguchi orthogonal designs offer minimizes the effect of aliases and measures error with minimum testing. In this context, aL9 orthogonal array of 3 factors and 3 levels was chosen (see Table 4).

The present study used ANOVA to determine the optimum combination of process parameters more accurately by investigating the relative importance of each parameter [29] . Table 5 presents the results of ANOVA for the tooth hardened depth (D_{t}). It is observed from the results (Table 5 and Figure 8) that the scanning speed (68.66%) is the most significant parameter, followed by power (25.14%). The rotary speed has the smallest effect (3.02%) in hardened depth. Statistically, the F-test determines whether the parameters are significantly different. A larger F value shows a greater impact on the machining performance characteristics [29] . Larger F- values are observed for scanning speed, as 25.24, and for power, as 7.18.

As seen from the ANOVA results in Table 6 and Figure 9, the influence of the scanning speed (66.06%) in the hardened depth of the flank is significantly larger. The power (27.28%) is the second most significant factor. Again, the rotary speed has the least effect (3.45%) on D_{f}. It is also observed that there is an error contribution of 3.21% in the hardened depth on the flank.

7. Conclusion

In this paper, an integrated approach used to build a hardness profile prediction model for AISI 4340 spline shafts heat treated by laser is presented. Numerical simulation carried out through 3D finite element model using Comsol Multiphysics software is discussed. A commercial 3 kW Nd:Yag laser system, a structured experimental

Figure 8. Effect of the parameters on the depth (D_{t}).

Figure 9. Effect of the parameters on the depth (D_{f}).

Table 4. Experimental planning.

Table 5. Variance analysis case of depth (D_{t}).

Table 6. Variance analysis case of depth (D_{f}).

design and confirmed statistical analysis tools are used to conduct the experimental study for the prediction model calibration and validation. The results reveal that the numerical simulation can effectively lead to a consistent and accurate model and provide an appropriate prediction of the hardness profile attributes under variable hardening parameters and conditions. With a maximumerror less than 15%, the validation process shows great concordance between predicted and experimental results.

Cite this paper

Mahdi Hadhri,Abderazzak El Ouafi,Noureddine Barka, (2016) Hardness Profile Prediction for a 4340 Steel Spline Shaft Heat Treated by Laser Using a 3D Modeling and Experimental Validation. *Journal of Materials Science and Chemical Engineering*,**04**,9-19. doi: 10.4236/msce.2016.44002

References

- 1. Kannatey-Asibu Jr., E. (2009) Principles of Laser Materials Processing. Vol. 4, John Wiley & Sons, Hoboken, 860.

http://dx.doi.org/10.1002/9780470459300 - 2. Steen, W., Watkins, K.G. and Mazumder, J. (2010) Laser Material Processing. Springer, London.

http://dx.doi.org/10.1007/978-1-84996-062-5 - 3. Chattopadhyay, R. (2004) Laser Assisted Surface Engineering Processes. Advanced Thermally Assisted Surface Engineering Processes, Kluwer Academic Publishers, Dordrecht.
- 4. Goia, F. and de Lima, M. (2011) Surface Hardening of an AISI D6 Cold Work Steel Using a Fiber Laser. Journal of ASTM International, 8, 315-318.

http://dx.doi.org/10.1520/JAI103210 - 5. Tobar, M.J., álvarez, C., Amado, J.M., Ramil, A., Saavedra, E. and Yáñez, A. (2006) Laser Transformation Hardening of a Tool Steel: Simulation-Based Parameter Optimization and Experimental Results. Surface and Coatings Technology, 200, 6362-6367.

http://dx.doi.org/10.1016/j.surfcoat.2005.11.067 - 6. Komanduri, R. and Hou, Z.B. (2001) Thermal Analysis of the Laser Surface Transformation Hardening Process. International Journal of Heat and Mass Transfer, 44, 2845-2862.

http://dx.doi.org/10.1016/S0017-9310(00)00316-1 - 7. Badkar, D.S., Pandey, K.S. and Buvanashekaran, G. (2011) Parameter Optimization of Laser Transformation Hardening by Using Taguchi Method and Utility Concept. The International Journal of Advanced Manufacturing Technology, 52, 1067-1077.

http://dx.doi.org/10.1007/s00170-010-2787-z - 8. Grum, J. and Kek, T. (2004) The Influence of Different Conditions of Laser-Beam Interaction in Laser Surface Hardening of Steels. Thin Solid Films, 453, 94-99.

http://dx.doi.org/10.1016/j.tsf.2003.11.177 - 9. Kar, A. and Mazumder, J. (1989) Three-Dimensional Transient Thermal Analysis for Laser Chemical Vapor Deposition on Uniformly Moving Finite Slabs. Journal of Applied Physics, 65, 2923-2934.

http://dx.doi.org/10.1063/1.342739 - 10. Ion, J. (2005) Laser Processing of Engineering Materials: Principles, Procedure and Industrial Application. Butterworth-Heinemann, Elsevier Butterworth-Heinemann, Burlington.
- 11. Komanduri, R. and Hou, Z.B. (2004) Thermal Analysis of Laser Surface Transformation Hardening Optimization of Process Parameters. International Journal of Machine Tools and Manufacture, 44, 991-1008.

http://dx.doi.org/10.1016/j.ijmachtools.2004.01.011 - 12. Reséndiz-Flores, E.O. and Saucedo-Zendejo, F.R. (2015) Two-Dimensional Numerical Simulation of Heat Transfer with Moving Heat Source in Welding Using the Finite Pointset Method. International Journal of Heat and Mass Transfer, 90, 239-245.

http://dx.doi.org/10.1016/j.ijheatmasstransfer.2015.06.023 - 13. Yang, J., Sun, S., Brandt, M. and Yan, W. (2010) Experimental Investigation and 3D Finite Element Prediction of the Heat Affected Zone during Laser Assisted Machining of Ti6Al4V Alloy. Journal of Materials Processing Technology, 210, 2215-2222.

http://dx.doi.org/10.1016/j.jmatprotec.2010.08.007 - 14. Toyserkani, E., Khajepour, A. and Corbin, S. (2004) 3-D Finite Element Modeling of Laser Cladding by Powder Injection: Effects of Laser Pulse Shaping on the Process. Optics and Lasers in Engineering, 41, 849-867.

http://dx.doi.org/10.1016/S0143-8166(03)00063-0 - 15. Safdar, S., Li, L. and Sheikh, M.A. (2007) Numerical Analysis of the Effects of Non-Conventional Laser Beam Geometries during Laser Melting of Metallic Materials. Journal of Physics D: Applied Physics, 40, 593-603.

http://dx.doi.org/10.1088/0022-3727/40/2/039 - 16. Sheikh, M.A. and Li, L. (2010) Understanding the Effect of Non-Conventional Laser Beam Geometry on Material Processing by Finite-Element Modelling. Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, 224, 1061-1072.

http://dx.doi.org/10.1243/09544062jmes1745 - 17. Rozzi, J.C., Pfefferkorn, F.E., Incropera, F.P. and Shin, Y.C. (2000) Transient, Three-Dimensional Heat Transfer Model for the Laser Assisted Machining of Silicon Nitride: I. Comparison of Predictions with Measured Surface Temperature Histories. International Journal of Heat and Mass Transfer, 43, 1409-1424.

http://dx.doi.org/10.1016/S0017-9310(99)00217-3 - 18. Patwa, R. and Shin, Y.C. (2007) Predictive Modeling of Laser Hardening of AISI5150H Steels. International Journal of Machine Tools and Manufacture, 47, 307-320.

http://dx.doi.org/10.1016/j.ijmachtools.2006.03.016 - 19. Rozzi, J.C., Pfefferkorn, F.E., Incropera, F.P. and Shin, Y.C. (1998) Transient Thermal Response of a Rotating Cylindrical Silicon Nitride Workpiece Subjected to a Translating Laser Heat Source, Part I: Comparison of Surface Temperature Measurements with Theoretical Results. Journal of Heat Transfer, 120, 899-906.

http://dx.doi.org/10.1115/1.2825909 - 20. Skvarenina, S. and Shin, Y.C. (2006) Predictive Modeling and Experimental Results for Laser Hardening of AISI 1536 Steel with Complex Geometric Features by a High Power Diode Laser. Surface and Coatings Technology, 201, 2256-2269.

http://dx.doi.org/10.1016/j.surfcoat.2006.03.039 - 21. Orazi, L., Liverani, E., Ascari, A., Fortunato, A. and Tomesani, L. (2014) Laser Surface Hardening of Large Cylindrical Components Utilizing Ring Spot Geometry. CIRP Annals—Manufacturing Technology, 63, 233-236.

http://dx.doi.org/10.1016/j.cirp.2014.03.052 - 22. Benedict, G.F. and Eskildsen, J. (1985) Method and Apparatus for Laser Gear Hardening. US Patent No. 4,539,461.
- 23. Zhang, H., Shi, Y., Xu, C.Y. and Kutsuna, M. (2003) Surface Hardening of Gears by Laser Beam Processing. Surface Engineering, 19, 134-136.

http://dx.doi.org/10.1179/026708403225002595 - 24. Pretorius, T. and Vollertsen, F. (2009) Simulation of the Distortion Manipulation of Gear Wheel Teeth by Thermal Pre-Stressing. Materialwissenschaft und Werkstofftechnik, 40, 479-484.

http://dx.doi.org/10.1002/mawe.200900480 - 25. Majumdar, J.D. and Manna, I. (2015) Laser Surface Engineering. In: Nee, A.Y.C., Eds., Handbook of Manufacturing Engineering and Technology, Springer, London, 2639-2676.

http://dx.doi.org/10.1007/978-1-4471-4670-4_27 - 26. Kakhki, M.E., Kermanpur, A. and Golozar, M.A. (2009) Numerical Simulation of Continuous Cooling of a Low Alloy Steel to Predict Microstructure and Hardness. Modelling and Simulation in Materials Science and Engineering, 17, Article ID: 045007.

http://dx.doi.org/10.1088/0965-0393/17/4/045007 - 27. Carlone, P., Palazzo, G.S. and Pasquino, R. (2010) Finite Element Analysis of the Steel Quenching Process: Temperature Field and Solid-Solid Phase Change. Computers & Mathematics with Applications, 59, 585-594.

http://dx.doi.org/10.1016/j.camwa.2009.06.006 - 28. Yung, K.C. and Zhang, B. (2010) Analysis of Process Parameters of Laser Structuring with Taguchi Method. Applied Physics A, 101, 385-392.

http://dx.doi.org/10.1007/s00339-010-5832-8 - 29. Mason, R.L., Gunst, R.F. and Hess, J.L. (2003) Statistical Design and Analysis of Experiments: With Applications to Engineering and Science. John Wiley & Sons, Hoboken.

NOTES

^{*}Corresponding author.