On the Effect of Distributed Regenerative Delay on the Stability Lobe Diagrams of Milling Processes

Regenerative machine tool chatter is investigated for milling operations with helical tools. The stability of a two-degreesof-freedom milling model is analyzed, where the cutting-force is modeled as a force system distributed along the rake face of the tool. Introducing a distributed force system instead of a concentrated cutting-force results in an additional short, periodically varying distributed delay in the governing equations of the system. It is shown that the additional delay significantly affects the stability of the machining operation, especially at low spindle speeds. This phenomenon is referred to as the short regenerative effect, and is studied by computing the stability lobe diagrams of milling operations via the semi-discretization technique. The sensitivity of the stability charts to the shape of the force distribution and the contact length between the chip and tool is investigated.


Introduction
The occurrence of harmful vibrations (chatter) during metal cutting processes is an important problem in manufacturing technology.Machine tool chatter has many unfavorable effects: it reduces the productivity and the surface quality, causes noise, enhances tool wearing, and even leads to tool damage in some cases.Therefore, avoiding or suppressing chatter is highly important for the machine tool industry.In the past few decades, a significant amount of research was conducted to find the governing physical phenomena behind chatter in order to understand its nature and describe methods to avoid or suppress it.
One of the most widely accepted explanations of machine tool chatter is the theory of surface regeneration introduced by Tobias [1] and Tlusty [2]: the vibrating tool leaves a wavy surface behind, which results in a varying cutting-force in the consecutive cut.Accordingly, delay effects appear in the models of metal cutting operations since the cutting-force exciting the tool motion is determined by the chip thickness, which depends both on the actual tool position and the delayed position at the previous cut.Hence, machine tool vibrations can be described using delay-differential equations, and the regenerative machine tool chatter can be considered as the manifestation of self-excited oscillations in a time-delay system.
Following the works of Tobias and Tlusty, a large effort has been put lately into the more accurate modeling of machine tool chatter.In this paper we follow the model of [3], where the stability of turning processes is investigated taking the distribution of the cutting force along the tool's rake face into account.The concept was experimentally verified in [4].Extension to interrupted turning operation was presented in [5].Here, we will extend the distributed cutting-force model to milling processes and perform the stability analysis.
The outline of the paper is the following.Section 2 presents the mechanical model of the system, the description of the cutting-force expression, and the final form of the governing equation.In Section 3, the semi-discretization technique [6] is applied to analyze the stability of the system numerically.Some numerical issues regarding computational efficiency are also highlighted.Section 4 considers the special cases of the investigated milling model, and a comparison to previously published results is given.Finally, the main results of this paper are shown in the form of stability lobe diagrams in Section 0, where the sensitivity to the size of the contact region between chip and tool as well as to the shape of the cutting-force distribution is also investigated.

Mechanical model
In this work we analyze the stability of end milling operations with helical tools.The corresponding mechanical model can be seen in Fig. 1.The model construction follows [7], where single-and two-degrees-of-freedom milling models are discussed assuming concentrated cutting-force.In many machining operations, including milling, the most compliant part of the machining system is the machine tool.In this case, we consider the workpiece as rigid, whereas the mill itself is modeled as a cantilever beam.This allows us to describe the dynamics of the machine tool-workpiece system by a twodegrees-of-freedom model.Accordingly, the equation governing the motion of the tool assumes the form where M, C and K denote the 2 × 2 modal mass, damping and stiffness matrices, respectively, and are the vector of general coordinates and the cutting-force vector.
Let us formulate an expression for the force components F x (t) and F y (t).Since a helical tool is used, the cutting force varies in the axial direction z along the cutting edges.Therefore, the cutting-force can be calculated by dividing the mill into elementary disks with height dz and integrating the elementary cutting-force acting on them.Hence we can write where a p is the axial depth of cut.The quantities dF x (t, z) and dF y (t, z) denote the elementary cutting-force components acting at axial immersion z, and they can be obtained as the sum of forces exerted on the different teeth of the tool, where j = 1,2, ..., N is the index of the teeth assuming an N-fluted tool.
Let us model the cutting-force on each tooth as the resultant of a force system distributed along a finite contact region on the rake face, that is, along the interface between the chip and the cutting edge.The length of the contact surface on the jth flute at a given axial coordinate z is denoted by l j (t, z).Note that the contact length is time dependent.In milling operations the cutting is intermittent, since the teeth of the mill enter and exit the workpiece material repeatedly.Therefore, the contact length l j (t, z) varies periodically.As the flute engages in cutting, the contact region evolves gradually over a finite amount of time, and then it drops down to zero when the cutting edge leaves the material.
Since the particular cutting edge segments of the different teeth do not enter the material at the same time, the contact length depends on j and z.However, now we assume that apart from a shift in time, the contact length l j (t, z) evolves according to the same law regardless the tooth index j and the axial position z.We suppose that for every cutting edge segment it increases linearly from 0 to the final value l during σ time, and then it remains constant until the cutting edge exits the workpiece.The linear increase is valid if the chip slips along the rake face with constant speed.Given the tool diameter D, the angular position φ j (t, z) of the jth cutting edge at axial immersion z, the entrance and exit immersion angles φ en and φ ex , the corresponding contact length reads (2) where mod is the modulo function.We can calculate the location φ j (t, z) of an elementary cutting edge segment using the function where Ω denotes the spindle speed, which is traditionally measured in rpm.The function Ψ(z) is the twist of the particular cutting edge.Assuming a uniform helix angle β, we can express the twist in the form where l p = Dπ / (N tan β) is the helix pitch [8].The entrance and exit immersion angles depend on the direction of rotation of the tool.The case when the cutting edge moves opposite to the feed direction at the entrance position is called up-milling, for which where a e is the radial immersion, and a e / D is called radial immersion ratio.
Whereas in case of down-milling, when the tool counterrotates and the flute enters the material in the feed direction, the angles become The contact length variation (5) in case of up-milling is illustrated in Fig. 2.
Knowing the contact length, we choose the local coordinate s Î [−l j (t, z), 0] to describe the force distribution along the rake face, see Fig. 1.Denoting the distributed force components by P j,x (t, z, s) and P j,y (t, z, s) we can write In order to treat the magnitude and the shape of the cutting-force distribution separately, we use the multiplicative decomposition ,y The function W j (t, z, s) characterizes the shape of force distribution, and assumes a normalized form such that W t z s s j l t z j ( , , ) .
( , ) The expressions d T F t z s j x , ( , , ) and d T F t z s j, ( , , ) y account for the cutting-force magnitude, and will be given later in this section.Prior to that, for convenience, let us transform the spatial description of the cutting-force into time.Assuming a constant spindle speed Ω, the cutting speed v = DπΩ/60 is also constant, and it is reasonable to assume that the chip slips along the face of the mill teeth with the same speed.Hence let us introduce the temporal coordinate θ = s/ν, which characterizes the time that elapses during a certain particle of the chip slip along the active face of the tool.The domain of θ is θ Î [−σ j (t, z), 0], whence σ j (t, z) = l j (t, z)/ν.According to Eq. ( 5), we can write Recall that σ denotes the time it takes for a particle of the chip to slip along the tool face.Using the temporal coordinate θ, we can rewrite Eqs. ( 10)- (11) in the following form where w j (t, z, θ) = νW j (t, z, νθ).According to Eq. ( 12), the following criterion applies to w j (t, z, θ): , )d .
Now we need to choose the function w j (t, z, θ), which characterizes the shape of cutting-force distribution along the rake face.
We can make a proper choice on w j (t, z, θ) via taking experimental data into account.In the past decades several measurements were performed to determine the distribution of stresses along the rake face of the tool.According to [9][10][11][12][13], there are two kinds of widely accepted shapes for both the normal and the shear stress distributions.Namely, some of the above measurements showed that the normal stress peaks at the tool tip and decays exponentially to zero at the end of contact.Alternatively, the other experiments reported a plateau in the normal stress distribution near the tool tip, and an exponential decay to zero after a certain length of the plateau.Whereas the experimental data on the shear stress distribution either followed the plateau-and-decay tendency, or featured a function that starts from a small value at the tool tip, increases to a maximum, and then decays to zero.In our model, we will approximate the latter function by a half sine wave, and we will also analyze exponential and constant distribution shapes to investigate the exponential and the plateau-and-decay distributions.In addition, considering Dirac delta force distribution yields the special case of concentrated cutting-force, which was analyzed in [7].Therefore, in this paper we assume that w j (t, z, θ) takes one of the following forms ( , Using these functions we assume that the shape of the force distribution is essentially the same for the different cutting teeth and axial immersions.The dependency of w j (t, z, θ) on j, t and z reflects only the fact that the same shape is scaled according to σ j (t, z).Therefore, we can introduce the scaling θ θ  = σ j t z ( , ), , and w tz t z, t z j j j θ σ ω σ θ ( ) = ( , ) ( , ( , ) ) for the cases where σ j (t, z) ≠ 0. This way we can represent the pure shape of the force distribution without considering the contact length.From (16) we get Besides, let us define w θ ( ) in a way that w θ ( ) ≡ 0 for 0 , which implies that the cutting-force is zero outside the contact region.Functions ( 17) can be seen in Fig. 3. Now let us focus on the force magnitude d T F t z j x , ( , , ) θ and d T F t z j, ( , , ) y θ .In machining the cutting-force magnitude is influ- enced by several factors, among which the most important ones are technological parameters like the depth of cut, the cutting speed, and the chip thickness.In the last century many cutting experiments were performed with the purpose of exploring the relation between the cutting-force magnitude and the above technological parameters.
Fig. 3 The investigated force distribution shapes along the rake face.
The most well-known expressions are the power law introduced by Taylor [14] and the cubic polynomial used by Shi & Tobias [15], but other kinds of force characteristics also exist, see e.g.
where subscripts t and r indicate the tangential and radial components, respectively.We can calculate the tangential and radial cutting-force components based on the above mentioned empirical relations with the technological parameters.In this work we will use Taylor's formula, the superscript T stands for Taylor as well.The Taylor expression reads where K t and K r are the tangential and radial cutting-force parameters obtained via measurements, and h j (t, z, θ) is the chip thickness along the rake face of a particular cutting edge.The cutting-force exponent q is usually assumed to be q = 3/4 leading to the well-known three-quarter rule.Note that sometimes a factor of ν −0.1 is also included on the right-hand side of formulae (19) in order to account for the dependence on the cutting speed.However, we neglect this term in our model.The coefficient g j (t, z, θ) is a screen function, which determines whether the jth cutting edge is engaged in cutting at axial immersion z: ( ( , )mod , Now the expression of all important quantities are available, and finally we can calculate the cutting-force F(t).Substituting ( 4), ( 14), (18), and ( 19) successively back into (3) yields ) Finally, let us formulate the chip thickness as a function of the tool position.According to the theory of regenerative machine tool chatter, as the cutting edge vibrates, it leaves a wavy surface behind.In the next revolution, the waviness of the surface changes the chip thickness, and hence modifies also the cuttingforce, which excites the tool motion.This results in the regeneration of the surface waves and the emergence of self-excited vibrations in the tool-workpiece system, also known as the machine tool chatter.Therefore, as the tool vibrates relatively to the workpiece, the actual chip thickness depends on the actual tool position and the position at the previous cut.Considering a circular tooth path approximation, the instantaneous chip thickness along the rake face of the jth cutting tooth reads It is important to remark that this formula neglects the displacement of the chip in the z-direction.In (22), f z denotes the prescribed feed per tooth, from which the actual feed per tooth differs due to the relative vibrations between tool and workpiece.
In fact, the parameter τ is called the regenerative delay, which now equals the tooth passing period τ = 60 / (NΩ).Correspondingly, σ is often referred to as short regenerative delay, and its influence on the system dynamics is called the short regenerative effect.In this work we assume that the ratio of the two delays is constant: σ = ετ.If we multiply this relation by the cutting speed ν = DπΩ / 60, it becomes clear that the constant ε is equivalent to the ratio of contact length l and the arc length Dπ / N between the neighboring teeth.Since this ratio is usually small, it explains the origin of the terminology short regenerative effect.
To sum up, the vibrations between the tool and the workpiece are governed by Eqs. ( 1), (21), and (22), which form a system of nonlinear delay-differential equations.The equations include a τ point delay, and an additional periodically varying distributed delay of maximum length σ.The kernel of the delay distribution is w j (t, z, θ), which originates from the shape of force distribution along the rake face of the tool.The coefficients of the differential equations are τ-periodic, whereas the period of the delay variation σ j (t, z) is Nτ.Thus the overall period of the system is τ.Let us assume that we can decompose the q(t) solution of the governing equations into a sum of a τ-periodic solution q p (t) = q p (t + τ) and a perturbation (t), The periodic solution q p (t) = (x p (t) y p (t)) T corresponds to the ideal chatter-free machining operation with constant feed per tooth f z , whereas the perturbation (t) = (ξ(t) η(t)) T represents the machine tool chatter.Substituting q(t) = q p (t) into the governing nonlinear differential equation yields the delay-free system where Q(t) is the τ-periodic function Equation ( 24) is an ordinary differential equation with τ-periodic forcing, hence it indeed has a τ-periodic solution, justifying that it is reasonable to assume the form (23). ( (25) Expanding F(t) into Taylor series with respect to q(t) around q p (t) and neglecting the higher-order terms yields the linearized system in the form The 2 × 2 matrix G j (t) is called the directional factor matrix given by where the components read , ( , ) cos ( , ) sin ( , ) ( , ) sin ,

sin ( , )cos ( , ). ϕ ϕ 3 Semi-discretization
In this section we use the first-order semi-discretization technique introduced in [6] to perform the stability analysis of Eq. (26) numerically.Note that other techniques can also be used to determine stability of systems with distributed delay, see, e.g., [17][18][19].As a result, stability lobe diagrams (or stability charts) are created, which identify the chatter-free domains in the plane of the spindle speed Ω and the axial depth of cut a p .In order to depict the stability charts, we use a fixed grid in the investigated region of the parameter plane, and in each gridpoint we determine the stability of the system via semi-discretization.
The main point of the semi-discretization technique is to approximate the solution operator of the infinite-dimensional delayed system by a large but finite-dimensional matrix.This is achieved through a piecewise constant approximation of the delayed terms in the governing equation.This way the delaydifferential equation is approximated by a sequence of discrete maps, and the stability analysis reduces to the eigenvalue computation of a linear map given by the so-called monodromy matrix, which describes the evolution of the system along its principal period τ.
The steps of the semi-discretization are discussed in [7].In that work, systems of the following form are considered with 0 and I indicating the 2 × 2 zero and identity matrices, respectively.Therefore, in order to apply the semi-discretization technique, we need to approximate (30) to reach a similar form as in (29).The two main differences between the above equations are the integral over z and the periodic variation of the length σ j (t, z) of the delay distribution in (30).Hence as the first step of the discretization, we replace the z-integral by a numerical quadrature using the rectangular approximation

∆
where m is the number of elementary disks the mill is divided into, and the following notations are used σ j,l (t, z) = σ j (t, z l ), W j,l (t, θ) = W j (t, z l , θ).The parameters of the numerical quadrature read z l = (l − 1/2)Δz and Δz = a p / m.
Secondly, we use a piecewise constant approximation of the delay variation σ j,l (t).The semi-discretization technique itself also approximates the delayed terms by taking them piecewise constant along the intervals t Î [ih, (i + 1)h], i Î .Here h = τ / p denotes the time step of the semi-discretization chosen in a way that the period τ is resolved into p intervals.In order to remain consistent with the semi-discretization technique, we approximate the σ j,l (t) delay variation in the same manner, namely we consider it as constant along the same t Î [ih, (i + 1)h) interval.Hence where the approximating constant delay σ j,l,i is defined as σ σ . = + ∫ 1 1  d System (30) with the approximations (32)-( 33) is now suitable for the semi-discretization method.Regarding the application of the semi-discretization we will highlight some important numerical issues, by which the computational time of the numerical stability analysis can significantly be reduced.
First of all, it is important to clarify that during semi-discretization we fix the period resolution p and not the time step h.This way it is possible to calculate several approximate quantities independently of the delay τ and hence of the spindle speed Ω.Therefore, we do not need to calculate them individually in each gridpoint of the stability charts, and a notable amount of computational effort can be saved.Accordingly, we introduce the dimensionless time t t , use the delay ratio instead of the delay σ j,l (t), and approximate it with ε j, l, i along ( ) ) in the following manner .
Note that introducing t allows us to express ( ) Consequently, using the relations σ = ετ, τ = 60 / (NΩ), and (13), it can be shown that becomes also τ-independent, as well as the approximation ε j, l, i .In a similar manner, the coefficient G j,l (t) = G j (t, z l ) in W j,l (t, θ) also has a τ-independent piecewise constant approximation G j, l, i .
Besides, the semi-discretization requires the approximation of the kernel function w j,l (t, θ) − w j,l (t, τ + θ) in W j,l (t, θ).It can be done using the integral we can write the w j,l,k,i approximation of (37) over t Î [ih, (i + 1)h) in the form Therefore, not only the periodic delay σ j,l (t) and the periodic coefficient G j,l (t) can be treated independently of τ, but also the kernel w j,l (t, θ).As a final remark, we must mention that the τ-independency could be reached only due to the fact p was fixed instead of h, and owing to a special property of the σ j,l (t), G j,l (t), and w j,l (t, θ) functions.Namely, as τ varies these functions should only be rescaled in time.Note that the dependency on the axial depth of cut a p cannot be transformed out in a similar manner.
Another important issue in the time-efficient computation of the stability charts is the construction of the monodromy matrix.In the course of semi-discretization, the evolution of the system over its principal period τ is given by p subsequent linear maps taking the system from t = ih to t = (i + 1)h, i = 0,1,..., p−1.Therefore, the monodromy matrix is constructed by multiplying p matrices.Since they are special sparse matrices, the monodromy matrix can be obtained in 'one step' as discussed in [7].In fact, we can build up the monodromy matrix recursively row-by-row using row operations instead of multiplying whole matrices.This also saves a considerable amount of computational time.
Finally, it is more reasonable to strive towards the computation of stability boundaries rather than determining the stability in each point of a large grid.The stability of the system is determined by the eigenvalues μ s of the monodromy matrix, also called as the characteristic multipliers.The system is asymptotically stable provided that each characteristic multiplier lies within the unit circle of the complex plane.Therefore, the stability boundary is defined by max .
Consequently, the stability analysis reduces to the root finding problem (39), where the semi-discretization technique provides the μ s eigenvalues.In order to find all solutions of (39), one can use the multi-dimensional bisection method introduced in [20], which is available in a Matlab package.This is an intervalhalving method initiated on a rough grid and supplemented by a continuation algorithm.With a sufficiently fine initial mesh the method is able to find all the stability boundaries. .

Special cases
In order to verify the results of the stability analysis, we can compare the stability charts to those of the previously studied cases with known solutions.Substitution of the Dirac delta kernel in (17) leads to the special case of concentrated cuttingforce.Whereas assuming β = 0, i..e., Ψ (z) ≡ 0 gives back the case of straight fluted tools.Besides, considering a symmetric system with M = m x I, C = c x I, K = k x I, and substituting G j (t) = G j,xx (t)I yields a single-degree-of-freedom milling model.Finally, if we assume a symmetric system with zero helix angle, consider a ( ) = constant coefficient instead of the periodically varying one, and keep the delay also constant σ j (t) ≡ σ, we end up with the single-degree-of-freedom model of turning.Running the semi-discretization code with the above assumptions, we checked the following special cases.First, the single-degree-of-freedom milling model with straight fluted and helical tools and the two-degrees-of-freedom model of milling with straight fluted tool assuming concentrated cutting-force in each case, which was also investigated in [7].Furthermore, we also computed the stability charts of the single-degree-of-freedom turning model with distributed cutting-force, which were derived analytically in [3].All these results were subject to comparison, and the stability charts obtained via the semi-discretization of system (30) agreed well with those of the special cases.
An example is presented in Fig. 4, where the first column shows the analytical results of [3] for the stability lobe diagrams of the single-degree-of-freedom turning model, whereas their numerical counterpart obtained using the semi-discretization is presented in the second column.The different rows of the figure assume Dirac delta, cosinusodial and sinusoidal force distributions in this specific order.The Dirac delta and the sinusoidal kernels are given in (17), while the cosinusoidal reads ( ) ( ) . The stability charts were computed on a 400x200 grid with the same parameters as those of the analytical stability lobes in [3].As can be seen, the agreement between the stability boundaries in the presented parameter range is good.

Stability lobes of milling
Figure 5 shows a series of stability charts for up-milling processes with radial immersion ratio a e / D = 0.05 using a two-fluted helical tool (N = 2) with helix pitch l p = 25 mm.The cuttingforce exponent is q = 0.75, whereas the cutting-force coefficients are K t = 107 • 10 6 N / m 1.75 and K r = 40 • 10 6 N / m 1.75 .The nominal feed per tooth is f z = 0.1 mm, thus the linearized cutting-force coefficients are K qf N m z q t − = ⋅ The corresponding undamped natural frequency and damping ratio of vibrations in both the x and the y directions are f n = 907 Hz and ζ = 0.02, respectively.The different rows of Fig. 5 correspond to the exponential, constant and sinusoidal kernels in (17), while the different columns assume delay ratios ε = 0.05,0.1 and 0.2.The special case of concentrated cutting-force with Dirac delta kernel is also presented at the top of the figure.
In the top right corner an unusually large range of spindle speed and axial depth of cut is depicted in order to show the stability lobes at the highest spindle speeds as well.In the meantime, in the top left corner and in the 3 × 3 stability chart series only that parameter region is shown, which is feasible and usual from technological point of view.As for the discretization, we resolved the period τ by p = 200 points, the parameters m = 100, N h = 20, N k = 100 were used, and the stability charts were computed on a 600x300 grid.
The resulting stability lobes exhibit some interesting properties.First of all, in [8] and [21] it was shown for a singledegree-of-freedom model that the mill helix induces additional unstable islands in the stability charts of milling processes.Moreover, it was also discussed that the helix-induced islands at highest spindle speeds are separated by horizontal lines where the axial depth of cut is an integer multiple of the helix (40) pitch.According to the top right corner of Fig. 5, these unstable islands appear in case of the two-degrees-of-freedom milling model as well, although they look like open curves when other stability lobes intersect them.Besides, the islands are still separated by horizontal lines (indicated by dashed lines) at integer multiples of the helix pitch.The stability lobe diagrams at the top also show that helix-induced islands appear at lower revolution numbers, too.
The effect of the cutting force distribution can be observed in the 3 × 3 stability lobe diagram series.The most important phenomenon owing to the short regenerative effect is the upward shift of the stability lobes at low spindle speeds.Sometimes this phenomenon is explained by the so-called process damping effect, see e.g.[22][23][24].The process damping is often introduced to account for the fact that the tool flank might interfere with the wavy surface of the workpiece.This way an additional force is produced, which pushes the tool away from the workpiece, and thus works as an additional damper with damping inversely proportional to the spindle speed Ω.This results in the upward shift of the stability lobes at low cutting speeds.An alternative explanation of the process damping uses the dependence of the instantaneous chip thickness on the vibration velocity of the cutting tool [13].However, based on Fig. 5 the expansion of the low-speed stable region can also be explained by considering a distributed cutting-force along the rake face instead of a concentrated one.It is important to note that contrarily to the process damping effect, the short regenerative effect also modifies noticeably the stability lobes at high spindle speeds.Therefore on the one hand, as the lobes change, the helix-induced islands tend to get outside the stable region and hence disappear from the stability charts.On the other hand, the stability boundaries intersect at lower depth of cut values, which suppresses some high-speed stable peaks of the stability charts.Therefore, we can conclude that although the additional short distributed delay stabilizes the system at low spindle speeds, it may also destabilize the process at certain high-speed parameter regions.
As a series of stability charts is presented in Fig. 5, it also reflects the sensitivity to the shape of force distribution characterized by ( ) ŵ θ , as well as to the length of the chip-tool interface measured by ε.The first two columns of the figure clearly show that the parameter ε, or equivalently, the contact length has a much stronger effect on the stability lobes than the shape of the force distribution.The short regenerative effect becomes more significant for large contact regions between tool and chip (large ε values), and is more pronounced in the case of the constant and the sinusoidal force distribution shapes.However, comparing the second and third column of Fig. 5 we can see only very slight differences between the related stability charts.It reveals the fact that above a certain limit value ε lim , an additional increase in ε does not change the stability lobes.
It can be explained as follows.As we increase ε, we assume that the contact region between chip and tool becomes larger.However, our model accounts for the fact that the chip needs a certain amount of time to slip along the active face.If the chiptool interface is supposed to be too long and the radial immersion is small, the chip does not have time to slip along the entire contact length, as the cutting edge leaves the material before the particles of the chip would reach the end of the contact region.In other words, the time the cutting edge spends in the material is not enough for the contact region to fully evolve.The size of the contact region cannot exceed a limit length -the arclength between the points where the cutting edge enters and exits the material.A division by the arclength Dπ / N between two successive cutting teeth gives the following limit value for ε: , N which is determined by the radial immersion ratio a e / D and the number N of cutting teeth.In our particular example the radial immersion ratio is quite small, a e / D = 0.05, and the mill is only two-fluted, N = 2.It yields ε lim = 0.144, which means that the assumption ε = 0.2 implies a contact length longer than what could actually evolve.Also, the condition in the third row of ( 5) and ( 13) is never fulfilled for ε > ε lim , the nonzero constant part of the contact length function in Fig. 2 is missing in this case.Therefore, above a limit value ε lim further increase in ε affects neither the actual contact length nor the system delay, i.e., results in the same stability charts for any ε > ε lim .Of course, this effect becomes significant only in case of very low radial immersions.

Conclusions
We considered a two-degrees-of-freedom mechanical model to describe machine tool chatter in milling operations with helical tools.We modeled the cutting force as a force system distributed along the tool's rake face, and applied the theory of regenerative machine tool chatter to compute the stability lobe diagrams of milling processes.The distribution of the cutting force results in the so-called short regenerative effect.We showed that the short regenerative effect influences the stability lobes at low spindle speeds in a similar manner as the process damping effect.Therefore, it can serve as an alternative explanation of the upward shift of the low-speed stability boundaries.
However, contrary to models with process damping forces, the short regenerative effect may alter stability lobes significantly also at high spindle speeds.As the two explanations yield some interesting differences in the behavior of the resulting stability charts, studying the structure of stability lobes in the case of a distributed cutting-force system and investigating the relevance of the short regenerative effect compared to the process damping effect will be subject of future research.

Fig. 1
Fig. 1 Two-degrees-of-freedom mechanical modell of milling operations with helical tool.

Fig. 2
Fig.2The variation of the chip-tool contact length and the corresponding length of the delay distribution. Fig.

Fig. 5
Fig. 5 Stabilty lobe diagrams of the two-degrees-of-freedom mechanical model of milling with helical tool. .