Mathematical Problems in Engineering

Volume 2013 (2013), Article ID 932013, 11 pages

http://dx.doi.org/10.1155/2013/932013

Research Article

Stability Prediction of Milling Process with Variable Pitch Cutter

Gang Jin,1,2 Qichang Zhang,1,2 Shuying Hao,³ and Qizhi Xie1,2

¹School of Mechanical Engineering, Tianjin University, Tianjin 300072, China

²Tianjin Key Laboratory of Nonlinear Dynamics and Chaos Control, Tianjin 300072, China

³School of Mechanical Engineering, Tianjin University of Technology, Tianjin 300384, China

Received 4 November 2012; Revised 25 February 2013; Accepted 25 February 2013

Academic Editor: Alexander Pogromsky

Copyright © 2013 Gang Jin et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

The use of variable pitch cutter is a known means to increase the stable limit depth of cut by disrupting the regenerative effect. In this paper, an improved semidiscretization algorithm is presented to predict the stability lobes for variable pitch cutters. Modeling efforts develop a straightforward analytical integral force model that can cover any case of piecewise continuous cutting regions regarding the helix angle. The proposed approach has been verified with the comparisons with prior works, time domain simulations, and cutting tests. In addition, the method is also applied to examine the effect of the tool geometries on the stability trends for variable pitch milling. Some new phenomena for certain combinations of parameters are shown and explained.

1. Introduction

Machine tool chatter is a self-excited vibration between the cutting tool and the work piece that would affect adversely both the surface finish and the dimensional accuracy of the workpiece and limit the efficiency of many machining operations. Furthermore, it promotes tool wear and may cause tool fracture and damage to the machine tool itself. In order to avoid chatter in machining, optimisation of machining conditions based on modelling and simulation of chatter is needed.

Generally, the chatter was usually analyzed and modeled in the frequency domain or in the time domain. The approach used to obtain the stability lobe diagram (SLD) is the most established method for predicting and preventing regenerative chatter. In the 1990s, Altintas and Budak [1] developed a stability method, known as the zeroth-order approximation (ZOA) method, by using the frequency domain and averaging the time-varying periodic coefficients of the derived delay differential equations (DDEs). However, this method cannot accurately predict the stability lobes at low radial immersion milling. To improve this case, a multifrequency solution was suggested to predict the stability by Merdol and Altintas [2]. Later, Insperger and Stépán [3–5] proposed the semidiscretization (SD) method, which is based on the discretization of the DDEs and approximates their infinite dimensional phase space by a finite discrete map in time domain. Bayly et al. [6] proposed a temporal finite-element analysis for solving the DDEs written in the form of a state space model and show good computational efficiency in low radial immersion. More recently, Ding et al. [7] proposed the full-discretization method (FDM) based on the direct integration scheme and verified that this method has high computational efficiency without loss of any numerical precision.

Meanwhile, a number of different strategies have been used to increase the stability of the machining system. Alternative methods focus on disrupting the regenerative effect to enlarge the stable zone of the SLD by expanding the stability frontier. For milling problems, the mechanism of regenerative chatter can be broken by using milling cutters with variable pitch. In case of variable pitch cutters, the phase between two waves is not constant for all teeth, thus disturbing the regeneration mechanism. Such disturbance reduces the modulation in chip thickness and slows down vibrations, which consequently increasing the stability of cutting. The effectiveness of variable pitch cutters in suppressing chatter vibrations during milling was first demonstrated by Slavicek [8]. He used a rectilinear tool motion for the cutting teeth and applied the orthogonal stability theory for irregular tooth pitch. By assuming an alternating pitch variation, he obtained a stability limit expression as a function of the variation in the pitch. Altintas et al. [9] adapted the analytical milling stability model to the case of variable pitch cutters, which can be used to analyze the stability with variable pitch cutters. For a given process and a desired spindle speed they optimized the tooth pitch by plotting the axial depth of cut as a function of the first and second pitch angle, which led to a stabilization of the process. Budak [10, 11] used a similar method for designing the variable pitch to maximize the stability limit. Note that the pervious analytical method is only suitable for variable pitch tools at high radial immersion. Recently, Sims et al. [12] investigated the stability of variable pitch and variable helix end mills by using adapted and time-averaged version of the SD method. Wan et al. [13] developed a unified method for predicting the stability lobes of milling process with multiple delays that arise from cutter run out or variable uniform pitched angle based on an updated version of the SD method. Sellmeier and Denkena [14] combined Ackermann’s method to control systems with delay with the method of the piecewise constant subsystems to investigate the stability of an unequally pitched end mill. Later, based on the model in [12], Yusoff and Sims [15] optimized variable helix tools combining SD method with differential evolution. However, the methodology presented in [12, 13] required discretization along the helical flutes to estimate the cutting forces, the calculation of the cutting force using axial slices will achieve less accuracy than integral relationships along the edge of cutter unless the number of axial slices is unrealistic to go to infinity. Dombovari and Stepan [16] recently introduced a general SDM-based mechanical model to predict the linear stability of special cutters with optional continuous variation of the helix angle. This model is suitable for considering the weight distribution of system regeneration. The shape of the weight distribution function was presented and discussed for nonuniform and harmonically varied helix angle cases. To the authors’ knowledge, beside the works in [16], there is a little work done to detect the effect of the tool geometries on stability trend.

With this issue in mind, an alternative approach to detect chatter for variable pitch milling is presented based on a straightforward model which captured cutting force through integral relationship. More specifically, the goal is to investigate the stability trend caused by tool geometries. The structure of this paper is organized as follows. In Section 2, a straightforward analytical force model is described. In Section 3, an updated SD method is utilized to obtain the chatter-free lobes. In Section 4, a series of experimental tests for validation are introduced. In the final section, the investigation of the influence of tool geometry on chatter stability is described.

2. Mathematical Model

The dynamic model of 2 DOF milling process considered here is shown in Figure 1. A cutter with radius and unequally spaced teeth rotating at a constant velocity Ω (rad/sec) is assumed to have flexibility in two orthogonal directions . The mathematical equation of this system is given by

where the vector contains the displacement of the tool tip in the and directions and the matrices , , and are the modal mass, damping, and stiffness matrices, respectively. Given the helix of the tool, the cutting force can be written as

where is the axial depth of cut and and are the differential force corresponding to an infinitesimal element thickness of the th tooth in the and directions, which can be, respectively, expressed as

Figure 1: Schematic mechanical model of two degrees of freedom milling system.

The angular position of the th tooth is defined by

where is the helix angle and represents the pitch angle between teeth (Figure 2). The tangential and radial differential cutting force components can be expressed as

where and are the respective linearized cutting coefficients along the tangential and radial directions of the tool. denotes whether the th tooth is cutting. This function is given by

where and are the start and exit angles of the th cutter tooth, respectively. The dynamic chip thickness can be expressed as

where is the pitch period corresponding to the pitch angle shown as in Figure 2(c).

Figure 2: Schematic drawing of helical end mill and workpiece. (a) The differential cutting force in the radial and tangential directions. (b) and (c) are the unfolded cutter circumference for uniform and variable helix cutter, respectively.

Substituting (3), (5), (6), and (7) into (2) and considering from (3), the resulting expressions can be yield as

where , and .

According to [17, 18], is piecewise continuous in an entry, middle-of-cut, and exit region, as shown in Figures 3(a) and 3(b). Based on (8), this feature is simultaneously controlled by , and . As shown in the bottom graphs of Figure 3, if denoting the intervals , and , the intersection of and is the integral interval for the current cutting region. Let such intersection be , such that

(a)

(b)

(c)

(d)

Figure 3: Schematic diagram of the three cutting regions, where the analytical force expression can be obtained as (top) (a) case one or (b) case two. Green regions represent the portion of the workpiece. Diagonal black lines represent the cutting edge of the tool. The two examples in the bottom part are the current cutting regions.

Note that (9) is effective only if cutting occurs. Thus, the introducing function , which determines whether the tooth is in contact with workpiece, can be defined as

Substituting (9) and (10) into (8) yields

where expresses the matrix of the directional cutting force coefficients and can be written as

The typical characteristics of the proposed model of cutting force can be listed as follows.

(1)Since the cutting force is captured along the flute in the direction of rotation, it is no longer dependent on instantaneous axial depth of cut anymore, but the rotation angle of the flute at the tool tip.(2)Compared with the prior cutting force models which consider the piecewise continuity throughout three regions of cutting (e.g., the work of Patel et al. [17]), (11) presents a unified cutting force expression suitable for any possible cutting regions. Therefore, the cutting force calculation following different regions is unnecessary, thus benefiting the application and programming.(3)The cutting force model can also be applied to uniform pitch cutter milling.

Combining (11) with (1), the governing equation of milling process for variable pith cutters considered the helix angle can be obtained. In the following section, methods are proposed to solve such equation to determine the stability lobes.

3. Stability Analysis

The periodic milling dynamics with time varying, self-excitation, and delay terms expressed in (1) are organized by using Cauchy transformation and are rewritten as

where

The nonlinear problem has been transformed to a linear problem. Noted that for linear damped system with single delay, the SD method proposed by Insperger and Stépán [3–5] can solve it well to determine the stability limits. However, given the appearance of more than one system delay attributed to the unequally spaced pitch for variable pitch cutter, the method is no longer suitable for present systems. For this reason, a new algorithm based on the SD method is proposed to deal with this problem directly.

Firstly, note that certain variables for variable pitch cutters are periodic with each tool revolution, which is in contrast to regular tools that are usually considered as periodic with each tooth pass. Thus, the cutter spindle rotation period is divided into number of discrete time intervals. Introduce symbol to represent the th time interval. Here, means the th time node. Then, the length of the interval can be calculated by

The number of intervals related to the delay item can be approximately obtained by

where indicates the operation that rounds positive number towards zero. Then, in , (13) can be approximated as

where

where means . and are the weighted factors that are used to relate with the state values at the two extreme node of the delayed interval .

For the initial value , the solution to (13) can be formulated as

The substitution of and into (19) leads to

where the associated matrices are given by

Here, denotes the identity matrix. Let

where

Combing (20) and (22), one can be recast into a discrete map form as

where each matrix is given by

The horizontal position of the discrete input matrices and in (25) depends on the value of corresponding to and begins from the column of and , respectively.

Based on (24) and (25), the following mathematical expressions can be established by coupling the solutions of the successive time intervals in period :

where is the Floquet transition matrix that gives the connection between and . The dynamic stability of the approximated system is evaluated in terms of the eigenvalues of the transition matrix , which is given by

Note that the policy proposed by Insperger and Stépán [4] can be used to decrease the computational burden by cancelling the rows and columns with values of zero related to every and . This cancellation will yield the -dimensional state vector

instead of the original -dimensional vector defined by (24). Correspondingly, the dimensional size of square matrix and will also be . For more detailed information, please see the works of Insperger and Stépán [4, 5]. Meanwhile, and will become 4 × 2 sized matrices and begin from the column of and , respectively.

4. Model Verification

4.1. Comparison with Prior Works

To illustrate the performance of the proposed approach on uniform and variable pitch milling tools, firstly, the works of Altintas and Budak [1] and Altintas et al. [9] are considered.

In Figure 4 the results using the proposed approach are compared to Altintas’ analytical method (zero order method), using the same parameter as the experimental work in [9]. It can be seen that the results of the proposed method agrees closely with those of analytical method for the uniform pitch cutter milling. The two methods obtain similar results for the unequally pitched end mill, except the spindle speed section from 8000 to 9000 rpm, in which the proposed approach obtains a higher stability chart. The unanimous conclusion was also gained from the work of Sellmeier and Denkena [14]. They checked the results of the point A and B in Figure 4(b) with the time domain simulations and found that point A is stable and point B corresponds to chatter. Moreover, they highlighted that the prediction error can be attributed to the replacement of time-varying system matrices with time-averaged matrices. In the following, experiments are conducted to capture this behavior as well as the experimental validation of the proposed method.

(a)

(b)

Figure 4: Comparison of the predicted stability lobes: (a) uniform pitch cutter with and (b) variable pitch cutter with .

4.2. Experimental Verification for Uniform and Variable Pitch Cutter

Cutting tests were carried out on a DMC 75V linear five-axis high-speed machining center. The experimental setup is shown in Figure 5. Two carbide end mills with  mm, , and were used. Both tools feature the same geometry except the tooth pitch. The teeth of tool number I are equally spaced (), and the teeth of tool number II are unequally spaced (). The end mill was mounted on an EPB-Monobloc HSK-A63 tool holder. The workpiece material was a block of aluminum alloy 6061-T651 of size 120 × 100 × 70 (mm) clamped on the Kistler dynamometer 9257B, which was adopted to measure the cutting force and was confined to the worktable.

Figure 5: Experiment setup.

The cutting force coefficients were obtained through a set of slotting experiments at constant axial depth of cut and different feed rates by using tool I. The measured tangential and radial cut force coefficients are 835.4 and 244.5 Mpa, respectively. In addition, as known from the works of Powell [19], the cutting force for a tool with variable pitch can be accurately predicted using the cutting force coefficients obtained from a regular tool. Thus, the same cutting force coefficients were used in the stability prediction for either uniform or variable pitch cutters.

The modal parameters at the tool tip were obtained through impact modal testing by using the commercial modal software LMS, and the result is listed in Table 1. One mode exited in the direction, whereas two modes exited in the direction. However, for the modes in the direction, only the second mode was considered in the proposed simulation method because of its higher flexibility relative to the first mode.

Table 1: Tool dynamic parameters.

Direction
Mode
Frequency,   

(Hz)
Stiffness,   

(× 10⁶ Nm⁻¹)
Damping ratios, (%)

1
1010.313
4.1343
5.95

1
711.512
9.4405
4.00

2
1018.700
3.9154
5.63

The stability predicted charts for tools I and II are shown in Figure 6. The stability lobes predicted by the analytical methods are also plotted for comparison. As seen from Figure 6, the stability lobes have a good consistency except for a deviation as shown in graph b of Figure 6, which has been framed by a brown rectangle. Two cutting tests are carried out at  rpm,  mm (i.e., 25% radial immersion ratio),  mm, and  mm/tooth for down milling (shown as the marks A and B1 in top graphs of Figure 6). One is used the regular pitch cutter, and the stability chart indicates that the milling system would be unstable under this condition. The measured cutting force (graph A-3) and its Fourier spectrum (graph A-4) [20] also prove that chatter vibrations are severe and evident. The chatter occurs at a frequency of 1002.3 Hz, which is close to the frequency of the first mode in the direction and the second mode in the direction of the machine-tool structure. The vibration and cutting force obtained from time domain simulations gave further evidence of chatter as shown in graphs A-1 and A-2. However, when the test was carried out using tool II under the same cutting condition, the chatter vibration disappears (graphs from B1-1 to B1-4 of Figure 6). Note that the time domain simulation takes into account the helix angle of the cutter, the time varying character of the system, and teeth jumping out of contact.

Figure 6: Comparison of stability lobes, forces, and surface quality obtained by experiments and simulation for uniform and variable pitch cutter.

Next, the cutting condition B2 of  rpm,  mm is adopted to evaluate the deviation of stability charts occurred at the spindle speed between 8000 and 10000 rpm. It can be seen from Figure 6 that the stability lobes obtained using the analytical method indicated that the milling system would be unstable under this condition. However, an opposite result is gained by the proposed method. In order to check the two different results further, cutting experiment and time domain simulation are conducted simultaneously. As observed in graphs from B2-1 to B2-4 of Figure 6, this test proves the instability for this cutting condition as well as the time domain simulation.

The previous experiments indicate that

(1)the proposed method is effective in stability prediction for uniform or variable pitch cutter; (2)for variable pitch cutters, the analytical method does result in a loss of precision in some region of stability lobes for variable pitch cutters. However, this case may be avoided for the proposed method.

5. Stability Chart Trends of Applying Cutter with Different Tool Geometries

Tool geometries have deep influence on the stability prediction for end milling. At a prescribed cutting speed range, tool geometries can be optimized to maximize stability limits against chatter for variable pitch cutter. Here, the stability trends are investigated as a function of tooth pitch variation, radial immersion ratios, number of teeth, and helix angle. Two types of tooth pitch variation presented in [21, 22] are adopted. They can be expressed as follows:

alternating pitch variation (APV):

linear pitch variation (LPV):

To express the varying strength of pitch for a variable pitch cutter, considering (29) and (30), variation factor is define as

In the following simulation, the modal parameters and cutting force coefficients applied were from the works of Altintas et al. [9]. Note that authors assume that such parameters do not vary with the changes in tool geometries.

For alternating variable pitch cutter, Figure 7 gives the concerned illustration for four different , that is, 0, 0.222 and 0.333. It can be seen that with increasing , larger maximum stability limits can be obtained, especially in high spindle speed domain. Meanwhile, the location of stability lobes move to the high speed domain, and the spindle speed corresponding to the maximum stability limits seemingly move to the opposite direction. This means that the policy to obtain maximum stability limits at a prescribed cutting speed range through adjusting is feasible.

Figure 7: Stability charts for alternating pitch variation cutter with , , and different under 50% radial immersion: (). (). ().

For linear variable pitch cutter, it can be seen from Figure 8 that the enlargement of could result in a significant improvement for absolute stability limit, especially for low speeds domain. However, the effect of on the stability gain becomes weaker and weaker when it is achieved to a relatively large level. In addition, while the minimum stability limits increase almost by 100%, the maximum stability limits decrease severely. This finding may indicate that such variation increases the stability gain at the expense of maximum stability limits. Thus, LPV is only an effect method for high-absolute stability limit.

Figure 8: Stability charts for linear pitch variation cutter with , , and different under 50% radial immersion: (). (). ().

At the fixed and for linear pitch variation cutter, as the small number of teeth is adopted, the absolute stability limits and stability area decrease evidently as shown in Figure 9. This means that although more number of teeth means bigger material removal rate (MRR) in general, the bigger MRR still can be obtained by fewer teeth due to the obvious effect of number of teeth on stability prediction for variable pitch cutter. Meanwhile, it can be observed that the location of stability lobes also moves to the left of the stability lobes diagram following the increase of number of teeth.

Figure 9: Stability charts for linear pitch variation cutter with , , and different number of teeth under radial immersion: (). (). ().

Figure 10 shows that at the fixed and , the difference among stability lobes caused by helix angles is slight expect for spindle speed around 8000–9000 rpm, in which the maximum stability limits decrease, and this trend becomes severer and severer following the increasing . It is worth noting that this is a specific phenomenon to variable pitch cutter, while no difference occurs for regular cutter as helix angle changes in high radial immersion.

Figure 10: Stability charts for alternating pitch variation cutter with (), , and different helix angle under 50% radial immersion.

6. Conclusions

In this paper, an improved algorithm for predicting the stability lobes of milling process with variable pitch cutter is presented. First, a straightforward analytical force model considering three piecewise continuous regions of the cutting for either uniform pitch cutter or variable pitch cutter is derived by integral relationship. Second, based on the cutting force modal mentioned above, the system governing equation of the milling process is obtained and solved by an updated SD algorithm. Experimental tests and time domain simulations are conducted to demonstrate the validation of the proposed method. Finally, research attention is paid on the influence of the tool geometries, for example, alternating pitch variation, linear pitch variation, the number of teeth, and helix angles on the stability trends in high radial immersion for variable pitch cutter milling. The following conclusions are achieved.

(1)The proposed method is effective for predicting stability lobes of milling process with either uniform or variable pitch cutter.(2)Analytical method does result in a loss of precision in some region of stability lobes for variable pitch cutter. However, this case may be avoided for the proposed method.(3)APV is a preferable policy to modify the tool design to look for the highest stability limits at a prescribed cutting speed range. For LPV, it is an effective method to increase the absolute stability limits in milling process. (4)The increasing number of teeth would lead to evident stability limits decrease and the movement of stability lobes to the lower speed domain. Compared the uniform pitch cutter [21], the effect of number of teeth on stability prediction is more obvious for variable pitch cutter.(5)The helix angle has a great influence on the stability limits for variable pitch cutters. This result is quite different from that for uniform pitch cutter. This means that optimizing helix angle of variable pitch cutter to improve stability limits leads to significant results.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Grant no. 11072168), Tianjin Research Program of Application Foundation and Advanced Technology (11JCYBJC05800).

References

Y. Altintas and E. Budak, “Analytical prediction of stability lobes in milling,” CIRP Annals, vol. 44, no. 1, pp. 357–362, 1995.

S. D. Merdol and Y. Altintas, “Multi frequency solution of chatter stability for low immersion milling,” Journal of Manufacturing Science and Engineering, vol. 126, no. 3, pp. 459–466, 2004.

T. Insperger and G. Stépán, “Semi-discretization method for delayed systems,” International Journal for Numerical Methods in Engineering, vol. 55, no. 5, pp. 503–518, 2002.

T. Insperger and G. Stépán, “Updated semi-discretization method for periodic delay-differential equations with discrete delay,” International Journal for Numerical Methods in Engineering, vol. 61, no. 1, pp. 117–141, 2004.

T. Insperger and G. Stépán, Semi-Discretization for Time-Delay Systems, Springer, New York, NY, USA, 2011.

P. V. Bayly, J. E. Halley, B. P. Mann, and M. A. Davies, “Stability of interrupted cutting by temporal finite element analysis,” Journal of Manufacturing Science and Engineering, vol. 125, no. 2, pp. 220–225, 2003.

Y. Ding, L. Zhu, X. Zhang, and H. Ding, “A full-discretization method for prediction of milling stability,” International Journal of Machine Tools and Manufacture, vol. 50, no. 5, pp. 502–509, 2010.

J. Slavicek, “The effect of irregular tooth pitch on stability of milling,” in Proceedings of the 6th International MTDR Conference, pp. 15–22, Pergamon Press, London, UK, 1965.

Y. Altintas, S. Engin, and E. Budak, “Analytical stability prediction and design of variable pitch cutters,” Journal of Manufacturing Science and Engineering, vol. 121, no. 2, pp. 173–179, 1999.

E. Budak, “An analytical design method for milling cutters with nonconstant pitch to increase stability, part I: theory,” Journal of Manufacturing Science and Engineering, vol. 125, no. 1, pp. 29–34, 2003.

E. Budak, “An analytical design method for milling cutters with nonconstant pitch to increase stability, part 2: application,” Journal of Manufacturing Science and Engineering, vol. 125, no. 1, pp. 35–38, 2003.

N. D. Sims, B. Mann, and S. Huyanan, “Analytical prediction of chatter stability for variable pitch and variable helix milling tools,” Journal of Sound and Vibration, vol. 317, no. 3–5, pp. 664–686, 2008.

M. Wan, W. H. Zhang, J. W. Dang, and Y. Yang, “A unified stability prediction method for milling process with multiple delays,” International Journal of Machine Tools and Manufacture, vol. 50, no. 1, pp. 29–41, 2010.

V. Sellmeier and B. Denkena, “Stable islands in the stability chart of milling processes due to unequal tooth pitch,” International Journal of Machine Tools and Manufacture, vol. 51, no. 2, pp. 152–164, 2011.

A. R. Yusoff and N. D. Sims, “Optimisation of variable helix tool geometry for regenerative chatter mitigation,” International Journal of Machine Tools and Manufacture, vol. 51, no. 2, pp. 133–141, 2011.

Z. Dombovari and G. Stepan, “The effect of helix angle variation on milling stability,” Journal of Manufacturing Science and Engineering, vol. 134, no. 5, Article ID 051015, 6 pages, 2012.

B. R. Patel, B. P. Mann, and K. A. Young, “Uncharted islands of chatter instability in milling,” International Journal of Machine Tools and Manufacture, vol. 48, no. 1, pp. 124–134, 2008.

B. P. Mann, B. T. Edes, S. J. Easley, K. A. Young, and K. Ma, “Chatter vibration and surface location error prediction for helical end mills,” International Journal of Machine Tools and Manufacture, vol. 48, no. 3-4, pp. 350–361, 2008.

K. B. Powell, Cutting performance and stability of helical endmills with variable pitch [Ph.D. thesis], University of Florida, Gainesville, Fla, USA, 2008.

T. Insperger, G. Stépán, P. V. Bayly, and B. P. Mann, “Multiple chatter frequencies in milling processes,” Journal of Sound and Vibration, vol. 262, no. 2, pp. 333–345, 2003.

E. Budak and Y. Altintaş, “Analytical prediction of chatter stability in milling, part I: general formulation,” Transactions of the ASME Journal of Dynamic Systems, Measurement, and Control, vol. 120, no. 1, pp. 22–30, 1998.

E. Budak and Y. Altintaş, “Analytical prediction of chatter stability in milling, part II: application of the general formulation to common milling systems,” Transactions of the ASME Journal of Dynamic Systems, Measurement, and Control, vol. 120, no. 1, pp. 31–37, 1998.