An analytical formulation of ZOA-based machining stability for complex tool geometries: application to 5-axis ball-end milling

The paper describes a computationally convenient analytical formulation of the stability of the cutting process with respect to self-excited vibrations in the case of five-axis milling based on the commonly used zero-order approximation. In the case of five-axis milling with general milling cutters, it is difficult to calculate stable machining process conditions for two main reasons. The first reason is the difficulty of calculating the mean value of the cutting force Jacobian with respect to the regenerative displacement (closely related to a milling directional matrix) for a generally inclined tool, and the second reason is the nonlinearity of this Jacobian with respect to the process parameters, which means that the problem cannot be reduced to a linear eigenvalue problem as is usual for linear cases (e.g. cylindrical milling with respect to the axial depth of cut). In the first part, this paper presents a modification of the calculation of the machining stability limits for a nonlinearly dependent cutting force Jacobian. A new formulation of this Jacobian for a general tool based on the surface integral over the tool and workpiece engagement region is presented which leads directly to the mean value of the Jacobian of the cutting force (direction matrix) without the need to calculate it as a function of time and then calculate the mean value over one revolution. The advantage is that if we can analytically describe the engagement area, we also obtain an analytical relation for the cutting force Jacobian. This is presented with the practical example of a generally inclined ball-end mill. This analytical formulation of the force Jacobian allows the calculation of its derivatives with respect to the technological parameters (depth of cut, step over, tilt and lead angles), which is useful both for the calculation of stability diagrams and for solving optimization problems related to machining stability.


Introduction
One of the main limits to efficient machining is tool and workpiece vibration, which can have various origins -from drives, their control [1] and the dynamic behaviour of moving machine parts [2] to the cutting process. Of the latter source, the most important is the phenomenon of self-excited vibration -chatter, which leads to rough surfaces and possible damage to the workpiece and tool, or at least to high tool wear. Chatter is a result of self-sustained oscillation due to tool workpiece interaction. The positive feedback that leads to the oscillation is caused by the fact that process stiffness depends on chip thickness regeneration, an effect firstly considered some 70 years ago by Tlusty and Polacek [3] for turning operation. In the case of milling, the problem is more difficult due to the time-dependent force response of the cutting process. Budak and Altintaş [4] developed an approximate method to solve chatter stability for cylindrical end-mills called zero-order approximation (ZOA). The method is based on averaging over one revolution of the time-dependent cutting force Jacobian matrix with respect to relative position change within two consecutive cuts, i.e. using only its zeroorder Fourier series term. It has been compared with more precise time-domain simulations and proven to be rather accurate for tools with standard geometry (tools without serration or variable helix angle) and operation (most importantly constant spindle speed) [5].
Milling stability calculation methods were primarily developed for cylindrical tools because of the simpler geometry of the process. There are two crucial differences between cylindrical and 5-ax ball-end milling which complicate stability calculation. For cylindrical tools, the cutting force Jacobian matrix is directly proportional to depth of cut (the directional matrix does not change with the depth of cut) and the process can be modelled in 2D and described by simple trigonometric formulas. In contrast, the cutting force Jacobian matrix for ball-end milling is not linearly dependent on the depth of cut and is difficult to calculate in 5-ax milling due to the difficulty of calculating the toolworkpiece engagement for non-zero tilt and lead angles. The first issue was addressed for 3-ax ball-end milling by Altintas and Lee [6]. Their approach was based on using an iterative calculation of the mean cutting coefficients with respect to variable depth of cut. The structural dynamics model of the tool and workpiece is considered in feed and normal directions without considering variable dynamics with the change of position.
Ozturk and Budak [7] and Ozturk et al. [8] proposed another method in their paper that successfully dealt with depth of cut nonlinearity. They also consider the tool's nonzero tilt and lead angles. The method is based on discretization of the engagement surface into layers perpendicular to surface normal direction and use of the layer thickness as a linear stability parameter. This calculation is done iteratively to obtain the correct limit depth of cut.
There are many approaches based on trigonometric relationships which work well in cylindrical mills, but due to the geometric complexity of the problem, the calculations or cutting force models are too specific to be useful in general 5-ax milling stability calculations, e.g. Wojciechowski et al. [9], Graham et al. [10], Hao et al. [11], Tsai and Liao [12], Shtehin et al. [13] or Lazoglu and Liang [14]. For this reason, it is often advantageous to use simulation software to calculate tool-workpiece engagement, as e.g. [15] or [16]. Chao and Altintas [17] used software determination of the engagement region in their article on chatter-free orientations. A similar approach was used by Li et al. [18] who followed the previous article and applied screw theory to define general tool-workpiece kinematics. Ju et al. [19] solid-analytical-based method, which Zhan et al. [20] used to analyse the effect of variable pitch, also uses software solutions for tool-workpiece intersection calculation.
However, the need to use a virtual machining tool to define the engagement slows stability lobe diagram (SLD) calculation or process optimization significantly. This article presents a method to efficiently (analytically) calculate ZOA-based stability lobe diagrams for ball-end milling, which offers an advantage because if contact between the cylindrical shank and workpiece is avoided, the shape of the contact is not affected by a change of tilt or lead angles.
In the analysis, the tool tip engagement situation (with non-zero undeformed chip thickness) is not considered due to the high level of cutting force model uncertainty near the tool tip linked to near-zero cutting velocity. A process damping model would be necessary for correct estimation of stability, particularly for vibration in the direction of the tool axis.
The first part presents a novel formulation of averaged cutting force Jacobian matrix with respect to relative toolworkpiece displacement based on surface integration over the engagement region. The surface integral is independent of the coordinate system so it can be parametrized in the workpiece coordinate system where the integral limits are simplest and as the engagement region is typically significantly smaller than the area covered by the rotating cutting edge (ca. half-sphere) it requires significantly fewer calculation steps than averaging time-dependent Jacobian. Furthermore, if the engagement surface can be described by (piecewise) analytical functions, the formulation of the Jacobian is also analytical, which enables its differentiation. This is important because it also allows to determine the derivatives of the stability exponent with respect to the technological parameters, which can be used for numerical calculation of SLD, distinguishing stable and unstable regions in the diagram or for optimization problems. Knowing the derivatives, it can be used, for example, to find the optimal orientation of the workpiece with respect to the process stability more quickly. The proposed approach can be considered as an analytical formulation and generalization with respect to the tool geometry of the approach proposed by Ozturk and Budak [7] and Ozturk et al. [8] for ball-end tools. The resulting Jacobian is generally nonlinear with respect to the process parameters (depth of cut, stepover) and it is therefore necessary to generalize Tlusty's prescription for these cases as well. This is the subject of the next section of the paper.
The second section presents a formulation of the nonlinear eigenvalue problem for constructing a stability lobe diagram for the limiting depth of cut (or other parameter, e.g. lead and tilt angles, step over) and the spindle speed, which generalizes Otto et al. [21] formulation of Tlusty's law for zero-order approximation (ZOA). This formulation is not limited to machining with a ball milling machine, but for example also to the side-milling process with a cylindrical tool, where instead of the axial depth of cut, with respect to which the force reaction is linear, we consider the radial depth, with respect to which the force reaction behaves nonlinearly. It is also shown how to calculate the derivatives of the parameters defining the stability limits given the analytical expression of the Jacobian of the cutting force.
In the third part, the new method is applied to ballend milling with general tool and workpiece orientation -the cutting process parameters are indicated in Fig. 1. The boundaries of this region are derived for the ball tool as analytical functions based on the engagement conditions expressed as geometric inequalities. This allows the mean value of the cutting force Jacobian with respect to the tool and workpiece regenerative displacement to be determined as an analytical function in the case of a ball milling machine. This allows derivations of this Jacobian with respect to the process parameters, which can be advantageously used to efficiently determine the SLD or to optimize the process with respect to the process with respect to the depth of cut, stepover or tilt and lead angles. A validation of the approach by comparing the separately averaged timedependent Jacobian cutting forces obtained numerically from a software simulation of material removal is performed at the end of the section. Although the method is derived from mathematical relationships, the validation ensures that there is no error in the numerical implementation.
The final sections present a numerical demonstration of the approach The approach is not restricted to depth of cut as a limit parameter, and therefore step over and lead and tilt angle stability lobe diagrams are also presented.

Novel analytical formulation of the cutting force Jacobian
In this section, an alternative formula for mean cutting force Jacobian is presented. The first part shows that generally, the averaging of integral based force modelling can be reformulated as a surface integral over the engagement area. The second part applies this to ball-end milling where the method leads to an analytical description of the force Jacobian needed for a stability lobe diagram calculation. The formula for the total cutting force is based on the commonly used Montgomery approach [22] of integration of cutting force differentials along the engaged cutting edge which are based on local process parameters. The force acting on an element of cutting edge of width dw is where [T] = [̂ ̂ ̂ ] is the transformation matrix between the local coordinate system at the cutting edge (LCS) and chosen global coordinate system (e.g. WCS) with a basis of the tangential ̂ (cutting velocity direction), normal ̂ (cutting edge created surface outer normal) and ̂ binormal vectors. The empirical model of the specific cutting force in the LCS is denoted as . For demonstration purposes, a shifted linear model identified on steel is used in the simulations below, i.e.
where e , c are edge and cutting force coefficients (represented as vectors in LCS basis). Its parameters (which depend on the cutting edge position) are the cutting velocity (magnitude) v c , inclination angle , rake angle and the dynamic undeformed chip thickness h which is the distance between a point on the cutting edge and the surface created in the previous cut. This distance is a result of the prescribed movement of the tool, which can be divided into translational and rotational movement of the tool as a rigid Fig. 1 Ball-end milling process geometry -step over a e , depth of cut a p and tilt and lead angles and coordinate systems used in the calcualtion: WCS -workpiece coordinate system, ECS -engagement coordiante system, LCSs -local coordinate systems at three diferent points at the cutting edge body relative to the workpiece, and dynamic disturbance, which is a result of deformations due to the compliance of the tool and the workpiece. The parameter can be approximated as where f t is feed per tooth related to tool centre point TCP , and F is a change of the tool axis orientation between two succesive cuts (see Fig. 2). This tool axis rotation can be either due to change of tilt and lead angle or change of the workpiece surface normal due to its curvature alont the toolpath. Vector denotes a point of the tool envelope which will be described more closely in the next paragraph. The symbol denotes the tool and workpiece regenerative dynamic displacement, see Fig. 2. It is calculated as a change of dynamic displacement of the tool centre point with respect to the workpiece position at the engagement within the time delay (given by time per one tooth revolution). The elemental chip width dw is the length of the cutting edge element projected onto a plane perpendicular to tangential direction.
In end milling, the force distribution on the workpiece and tool can be modelled as a point force which is obtained by line integration of all the force differentials expressed by Eq. (1) along all curves j representing j-th the tool cutting edge. It is assumed that all the curves lie on the same surface of revolution around the tool axis given by unit vector ̂Ω . This surface can be represented using Rodrigues' rotation formula with rotation matrix [R Ω ( )] corresponding to a rotation around Ω by an angle .
The LCS basis vectors on the engagement surface HGL(see Fig. 1) can be written in the following form

The transformation matrix [T] is composed of the basis vectors
The elementary chip width defined above as a projection of cutting edge element length to a direction perpendicular to the cutting velocity direction can be written using the LCS vectors as As pointed out by Otto et al. [21], the stability is based on linearization of the cutting force Jacobian with respect to regenerative dynamic displacement where g is a characteristic function which returns 1 if a point s on a j-th cutting edge is engaged and 0, if not.
The averaged cutting force Jacobian over one revolution necessary for the ZOA-based stability lobe diagram calculation is The previous formula for mean cutting force Jacobian contains double integration: one integral along the curve describing the cutting edge and a second over one revolution. The main point of this method of efficiently calculating the mean Jacobian is to transform the double integral in the specified variables into a general surface integral over an engagement surface. This requires division of the integrand by surface point distance from the tool axis of rotation ( ) , i.e.
=̂×̂. where the contact surface (engagement) between the tool and the workpiece is denoted S. The derivation of the formula is outlined in the Appendix A.1. The tool axis (axis of rotation) is defined by the vector ̂Ω . The cutting edges are space curves under the condition that the distance from the axis of rotation is constant at a given axial position An important observation is that although we started from a specific parameterization of the tool envelope, all inputs to the calculation are independent of this parameterization and work only with a general surface x, which can be parameterized arbitrarily to make the calculation as simple as possible. The functions in the integrand and parametrization of the engagement surface necessary for the mean force Jacobian calculation are defined in the following subsections. The most convenient parametrization of the tool-workpiece engagement surface is in the engagement coordinate system as used by Altintas and Lee [6] previously for zero lead and tilt angle milling.
In the case of the cylindrical milling tool geometry, it is possible to show the agreement of the formulation with the solutions in the literature. In the following section, an example of matrices for this case is given. The main advantage of this method is for more complex geometries where it is more difficult to define the belting of the tool. In the second part, the solution will be presented on a ball cutter when only the spherical part of the tool is belted. The implementation of a calculation including shank would require more complex definitions of the shanking conditions, the need to take into account possible different tool orientations in adjacent grooves, especially on curved surfaces, and the practical benefit to the design of the toolpaths would probably not be significant.

Application to cylindrical milling
The following section demonstrates the calculation on a simple case of a r diameter cylindrical milling cutter. For this case, the surface integral formulation does not provide any particular advantage, but due to the availability of solutions in the literature it is possible to compare the resulting matrices. Geometrically, this is a simple problem where the tool axis coincides with the normal to the workpiece surface. We can therefore work in the engagement coordinate system (ECS), which is given by feed F , cross-feed C and surface normal N directions, as shown in Fig. 3. In the first step, axis of rotation and parametrization of tool envelope surface is defined Using the formulas (4), the local coordinate system basis vectors with respect to the ECS are calculated: In order to calculate the surface integral, the integration limits must be determined in a given parametrization of the tool envelope. For explanatory purpose, this will be done in a manner used for the ball mill in the next section. The engagement surface points must satisfy the following conditions: (I) Point must be outside the volume removed in the previous row. The previous row is bounded by a plane, which is a moving cylinder envelope surface. This plane has normal vector C and its position is given This leads to integral limits st , ex for down-milling r cos ≤ a e − r,resp. for up-milling r cos ≥ −a e + r. (II) Envelope point must be inside the original workpiece volume also bounded by a plane Fig. 3(II). This plane has normal vector N and its position is given by the axial depth of cut a p .
This leads to 0 ≤ z ≤ a p (III) Condition of positive chip thickness on the tool envelope, i.e. ⋅ F > 0 , see Fig. 3(III). The condition leads to the following range of the angular parameter: 0 ≤ ≤ .
Taking into account the conditions (I) and (III), we have , ex = for down-milling and st = 0 , ex = arccos 1 − a e r for up-milling. The previous calculations of the integration limits could be computed based on a simple sketch, but the purpose of the example is to demonstrate the procedure of finding the integration limits using a system of inequalities, which will be used later in the more complex case of a ball mill.
The distance between tool envelope point and tool axis is for cylindrical mill equal its radius, ( ) = r . To allow comparison with the literature, the cutting coefficient vectors are expressed through tangential cutting coefficient K c and ratio of normal and tangential component K r in the form The result leads to the same Jacobian (directional matrix) as presented by Budak and Altintaş [4].

Tlusty's law for nonlinear ZOA-based cutting force Jacobian
Calculation of the stability lobe diagram requires finding a spindle speed dependence of a particular technological parameter for which the machining process is stable. Both structure and cutting process react by force to small displacements or velocity perturbation. Linearization of these responses leads to a dynamical system suitable for stability analysis. The structural response (mass, stiffness and damping matrices) can be modelled by FEA or the system can be described by FRFs at the tool-workpiece contact.
The aforementioned zero-order approximation (ZOA) simplifies the problem by averaging the time-dependent cutting process force interaction (linearized) between the tool and workpiece. This allows for relatively simple formulation of the problem in the frequency domain. After separating the static and dynamic parts of the motion equations (Floquet approach for milling, see Insperger and Stépán [23]), the machining stability problem is in the Laplace domain mathematically posed by a homogeneous equation is the structural transfer function (matrix). It is not necessary to specify the coordinate system of the structural transfer function matrix or the cutting force Jacobian matrix at this point, but it is necessary to ensure that they are related to the same coordinate systems. One of the convenient coordinate systems is the engagement coordinate system determined by a a machined surface normal at the tool location on the workpiece and feed direction. The use of the transformation of transfer functions from the coordinate system on the workpiece and on the tool to a common coordinate system in five-axis milling is described in more detail by Li et al. [18]. The argument of the transfer function, the Laplace parameter, has a real part, denoted c , which denotes the exponential evolution of the solution (i.e. the stability of the system being described), and an imaginary part, which corresponds to the frequency of oscillation and is denoted c . Since we are interested in finding the boundaries of the unstable domains in the machining parameter space where the real part of the Laplace variable is zero, it is sufficient to describe the system via the frequency response function (FRF), i.e. transfer function [Φ(i c )] resticted to purely imaginary Laplace parameters (frequencies). The vector (i c ) is the Laplace transform of the eigenvector of the self-excited oscillation (at the stability boundary). [∇ Δ ] is the cutting force Jacobian with respect to the regenerative displacement averaged over the tool revolution. This Jacobian is closely related to so-called directional matrix (see e.g. Budak and Altintaş [4]), which is used to calculate stability in cases where the force depends linearly on the depth of cut. The Jacobian depends on various process parameters, of which the depth of cut a p is considered in the analysis, step over a e , lead angle L , and tilt angle T . The time delay is the time between successive cuts. This homogeneous system of equations has nontrivial solutions only for certain combinations of parameters. In the following, the parameter chosen to calculate the lobe diagram from the above parameters a p , a e , L and T is denoted by p and the other parameters are omitted from the formulas for brevity. The condition for the existence of a non-trivial solution for a system of homogeneous equations is as follows The lobe diagram (boundaries of the unstable parameter regions) is calculated as a set of parametric curves (Ω( ( c )), p( c )) with frequency c as its parameter. The spindle speed directly related to the time delay .
Otto et al. [21] presented an approach to problems where the force Jacobian depends linearly on the depth of cut (or chip width). In that case, the stability condition can be transformed to the linear eigenvalue problem and the real limit depth of cut and corresponding real delays (and hence the spindle speed) are determined from the real and imaginary part of the eigenvalue.
The generally non-linear Eq. (16) can be treated similarly to the linear case. The auxiliary eigenvalue problem where = 1 − e −( c +i ) −1 is solved together with stability limit condition c = 0 , which is transformed to more direct condition Re = − 1 2 . It means we have a set of non-linear equations for the stability parameters p, which are solved for each chatter frequency c The time delay and related spindle speed are calculated from the eigenvalue using the same formulas as in the linear case The generalization is very similar to the formula proposed by Otto et al. [21] for systems where process stiffness is directly proportional to the depth of cut but it cannot be transformed to eigenvalue problem as in the linear case. There are various methods how to solve the equations numerically using for instance the Regula Falsi method or Newton-Ralphson method. In the latter, the analyticity of the proposed force Jacobian calculation can be utilized, as the calculation of the eigenvalue derivative p is straightforward. Knowing the eigenvalue and corresponding eigenvector for a given value of p, one can express the derivative as The derivative can be used for numerical calculation of the lobe diagram (Newton-Ralphson method) or for optimization problems. Another application is distinguishing between stable and unstable domains in the stability lobe diagram. As the problem is nonlinear, several values of the limit parameter for a given chatter frequency may exist. It is necessary to apply a method that is able to find all relevant limit parameters; for a possible approach, see [24].
The aforementioned summary and generalization of the previous results hold for ball-end milling as well as other tool geometries.

Application to ball-end milling
The geometry description generally follows notation as summed up by Ozturk et al. [25] and Li et al. [18]. This geometry is slightly simplified: the curvature of the machined surface is not considered and the toolpaths are (at least approximately) parallel. This however does not significantly limit its application to machining stability calculation which is considered (usually implicitly) for quasi-stationary conditions [26], which would not be possible on highly curved surfaces. The following forumulas are derived assuming that only the spherical part of the tool is engaged. To consider also the cylindrical shank, it would be necessary to divide the area into two domains and to solve the problem here separately. This would add another source of nonlinearities to the SLD calculation.
The preceding calculations at the general level do not require the specification of a coordinate system. However, for the application of Eq. 8 in relation to a specific tool, it is advisable to introduce convenient coordinate systems in which it is easier to introduce a parametrization of the engagement area for the calculation of the surface integral and to find the integration limits. These points will be the subject of the following subsections. The Wolfram Mathematica implementation of the Jacobian calculation presented in the sections bellow is in the Appendix A.4 Coordinate systems The coordinate systems used in the analysis respect the notation used by Li et al. [18]. The workpiece coordinate system (WCS), engagement coordinate system (ECS) and local coordiante systems (LCS) on the cutting edge are shown in Fig. 1. The engagement coordinate system is given by the feed direction vector ̂F , outer normal vector to the desired surface ̂N , and cross-feed direction ̂C normal to both, such that ̂F,̂C,̂N form a right-handed coordinate system.
The advantage of the ball-end tool is that the contact geometry does not depend on the tilt and lead angles (if only the spherical part is immersed in the workpiece). These angles are defined as angles of rotation about ̂C (lead angle) and ̂F (tilt angle) axes as illustrated in Fig. 1.
The point on the tool-workpiece engagement surface is conveniently expressed in the ECS independently on tilt and lead angles where the spherical surface is parametrized by the azimuthal angle and polar angle as shown in Fig. 4b.
As the origin of the ECS is placed at the centre of the spherical envelope, the tool envelope normals at any given point on the contact surface have the same direction as position vectors , i.e.
where r is the tool radius.
The inclined tool axis in the ECS is given by where T is a tilt angle and L is a lead angle. The minus is due to the clockwise spindle rotation typical for most machine tools.

Engagement region
Depth of cut a p is a maximum distance between the plane defining workpiece surface before material removal and engagement surface points in the ̂N direction, (20) = r sin sin ̂F + r cos sin ̂C − r cos ̂N and step over is the distance between adjacent toolpaths and is denoted as a e , see Fig. 4. It can be positive or negative to signify the direction of the uncut material with respect to the cross-feed vector. The following calculation of the engagement region is limited to situations when only the spherical part of the tool is engaged. Though the mathematical formulation is stated in more general form, the conditions are not significantly different from [27]. The engagement surface is given by the points on the tool envelope (spherical surface), which must satisfy the following conditions (I) To be outside the volume of material removed on the adjacent toolpath -approximated by a cylindrical surface (representing tool motion envelope) with an axis at a distance equal to the step over a e from the current toolpath, see Fig. 4I. The boundaries of the engagement surface are shown in Fig. 4.
The following calculations are done in order to express the above conditions in the parameters of the spherical engagement surface using analytic geometry. The boundary curve I (red) resulting from the first condition can be formulated as the intersection of the spherical surface (tool envelope) given by Eq. (21) and the cylindrical surface (volume removed by a spherical tool on the previous toolpath) defined parametrically as The following vector equation for intersection of a sphere (tool envelope) and a cylinder (adjacent toolpath envelope) must hold on the engagement boundary curve (red curve in Fig. 4)

This leads to the following limits in and
The condition II (semi-transparent yellow domain) formally defining the depth of cut can be formulated as which leads to the condition The last condition ⋅̂F ≥ 0 expresses the requirement of a positive chip thickness. In the given parametrization of the spherical surface, it leads to In summary, the conditions lead to the integration limits for the surface integral. These split into several cases depending on the condition in Eq.
the whole engagement surface is calculated as the sum of two surface integrals due to the condition in Eq. (23). The limits of the first integral represent boundaries of the half of spherical cap surface marked by the green semicircle in Fig. 4b. The limits are The second surface integral has the following limits The previous formulas hold, if the following is satisfied which means that there is not a gap between consecutive passes. If the condition does not hold, the process corresponds to slot milling with the limits (26) 10 = 0, One of the great advantages of the proposed approach is the analyticity of the formulation and hence the ability to calculate cutting force Jacobian derivatives with respect to process parameters. To calculate the derivatives with respect to step over and depth of cut, the derivatives of the limits need to be calculated.

Specific cutting force parameters
The specific cutting force coefficients are generally dependent on local cutting conditions which typically change along the cutting edge. It will be assumed that the cutting force coefficients depend on local rake angle (s) , inclination angle (s) and cutting velocity magnitude v c , where parameter s determines the position on the cutting edge.
To be able to incorporate the parametrically dependent cutting coefficient into the surface integral, it is necessary to be able to assign the parameter s to the parameters of the engagement surface , so that we have Let us assume that the cutting edge is parametrized by an azimuthal angle in the tool coordinate system. This leads to the following relationship between the parameters Due to tool rotation, the velocity magnitude v c can be obtained directly as the product of spindle speed Ω and the distance from the axis of rotation ( , ) The rake and inclination angles' dependence on parameter s is assumed to be an input based on tool parameters specified by its manufacturer.
The cutting force Jacobian and its derivatives The cutting force Jacobian under the previous assumptions on process geometry and specific force model is calculated from Eq. 8. In order to make the formulas more compact and readable, the following auxiliary function It may be physically interpreted as a stress gradient of the cutting process averaged over a tool revolution. The resulting formula for the Jacobian is , (s( , )), (s( , ))).  Although not explicitly stated due to the limited space of the text column, the auxiliary function is considered to be a function of the parameters s and , i.e.
[P] ≡ [P]( , s) , and the upper and lower limits are functions of , i.e.
2• ≡ 2• ( ) . This notation will be also applied in the following formulas, if not stated otherwise.
The significant advantage of the proposed approach is the ability to analytically calculate the derivatives of this Jacobian with respect to the process parameters, which speeds up both the calculations of the SLD and has applications for process optimization, where it allows to identify how which parameter changes the stability exponent through the process stiffness. Calculations of the derivatives of the Jacobian force with respect to the depth of cut and stepover, which only appear in the integration limits, are based on the Leibnitz integral rule The first formula is the contribution to the area integral from the distribution of the mean value of the force gradient per unit area along curve II in Fig. 4; the second formula is the contribution from the distribution along boundary I.
In order to calculate the derivatives of the Jacobian with respect to the lead and tilt angles for a spherical cutter, it is necessary to differentiate its integrand, i.e. the transformation matrix and the cutting coefficients (if the dependence on face angle, blade inclination or cutting speed is considered in the specific cutting force model).
The derivatives of c and derivatives of rake and inclination angle with respect to cutting edge parameter s come out of models of the tool geometry and specific cutting force. The last derivative comes from the cutting edge parameter definition, which might be for instance by Eq. (38). The formula for the derivative of the Jacobian with respect to the tilt angle can be formulated analogously. The derivative calculations are straightforward; however, due to their complexity, it is advisable to use symbolic mathematics software. An example of the calculation is presented in Appendix A.4. It is written with an emphasis on consistency with the method described above. For real applications, it is useful to perform e.g. analytical pre-calculation by normalizing vectors or expressing transformation matrices into functions separately. Calling multiple functions slows down the computation. For a rough idea of the speed of calculation, which may vary according to the computer used, we will give the calculation times we measured on our PC. The computation in Wolfram Mathematica presented in the appendix took about 0.7 s. After merging and simplifying the functions in Mathematica and transforming it to Matlab script, the computation of the Jacobian on the PC for different conditions takes about 0.03 s.

Numerical validation of the Jacobian calculation method
The above formulations of the Jacobian are based on mathematical identities and a mathematical description of the geometry of the process, and any objections to the validity of the resulting formulae should also be based on mathematical reasoning. However, given the complexity of the description, it is very easy to make mistakes (especially with respect to the orientations of the coordinate systems and the definitions of the technological parameters) and it is therefore advisable to validate the method with independent solutions on selected paths. The cutting force Jacobian calculated using the formula in Eq. 41 is compared with Jacobian calculated numerically using material removal software. The software is briefly described in Appendix A.2. The material removal that physically takes place on the individual teeth is replaced by a Boolean subtraction of the tool envelope from the workpiece represented by voxels in discrete positions m corresponding to the tool positions after each tool revolution (or after rotation by the tooth pitch). In Fig. 5a, there is a scheme of perturbed tool trajectory designed to numerically calculate the cutting force Jacobian. The dashed line represents unperturbed linear trajectory, black/grey discs unperturbed discrete tool positions and the coloured discs highlight perturbed tool positions in feed, cross-feed and surface normal directions. The actual feed vector calculated from the tool positions is split into the prescribed feed f t and perturbation according to Eq. 2. The position points distribution was chosen in a way that the prescribed feed is f t = 0.1 mm and the magnitude of the feed vector perturbation is d = ‖ ‖ = 0.01 mm . The numerical approximation of the Jacobian is based on cutting forces (m) calculated at the points m = {1, 2, 3, 4, 5, 6} as presented in Fig. 5a.
The cutting coefficients for the validation roughly correspond to orthogonal coefficients for carbon steel C45. The tool used for simulations (e.g. see slot milling in Fig. 5b) has a diameter 8 mm and 2 teeth. The unperturbed feed is chosen 0.1 mm and its perturbation is 0.01 mm. This leads to the perturbed feed (in Fig. 5c) and perturbed cutting forces shown in Fig. 5d-g, which are calculated according to methods presented in [28], for details see Appendix A.2.
The resulting time-dependent gradient is averaged and compared with gradients calculated by the proposed method. The comparison is in Table 1. The small difference can be attributed to the numerical error of the voxel discretization in the simulation and the fact that the toolpath perturbation affects not only the undeformed chip thickness but also the tool-workpiece engagement area.

Demonstration of the solution
A rather simple thin-walled plate dynamics simulated in FEA SW, where the transfer matrix was chosen to demonstrate the method. The FRF used for calculation is shown in Fig. 6 together with the dominant eigenshapes. The transfer matrix in ECS is available in the attachment as a Matlab file. The cutting coefficients are the same as in the previous section. The tool for validation has an 8-mm diameter and zero inclination (helix) angle.
The resulting lobe diagrams for various step over values, inclination, and lead angles are shown in Fig. 7. The points of the lobe diagram were plotted point-wise for each frequency without connecting the obtained solutions, giving a better idea of the solutions and their multiplicity. For clarity, areas of instability were additionally marked with grey fill.
Several observations can be made from the graphs. Firstly, it can be seen that the farther away the points of engagement are from the tool axis on average, the greater the critical depth of cut. This is consistent with the proposed reformulation of the integrand to calculate the mean cutting force Jacobian, where the distance from the axis of rotation is in the denominator (it can be interpreted that the smaller the distance on the engagement from the axis of rotation is in average, the longer the cutting edge spends in the engagement and higher cutting force reaction in average is). Secondly, the unusual lobe diagram for tilt 80 • clearly shows The physical explanation of this behaviour lies in the tangential and normal components of the cutting force. This can be seen in particular in Fig. 7f, g, where the same inclination in absolute value is chosen, but the opposite direction of the unmachined material with respect to the cross-feed direction, so that in the case of graph (f) the contributions of the tangential and normal force components are partially cancelled and in the case of graph (g) they add up, resulting in different machining stability. The effect can be also observed in cylindrical milling if a radial depth of cut is taken as its stability limit parameter. In the case of up-milling, the tangential and normal component projections into a direction perpendicular to feed direction act (on average) in the opposite direction and in the case of down-milling they act in the same direction. This effect has its analogy in cylindrical milling cutters, where due to the linearity of the Jacobian with respect to the axial depth of cut, this parameter is almost exclusively used for SLD drawing. The effect of radial depth of cut, with respect to which the Jacobian is nonlinear, or the effect of up or down milling is studied, for example, by Sanz-Calle et al. [29,30] through the stable axial depth of cut. The method presented here makes it possible to analyse directly the influence of parameters with respect to which the Jacobian of the cutting force is nonlinear. It was shown in Sect. 3 that almost any process parameter can be chosen as a limit parameter for machining stability. As a demonstration, stability lobe diagrams for spindle speed and lead angle, tilt angle and step over are shown in Fig. 8. SLD for the lead angle in Fig. 8a shows higher stability for the negative values of the lead angle. A small tilt angle was also chosen so that for small negative lead angles the tool tip does not get into the engagement, because there is no  Fig. 7, can be explained by the projection of the tangential and normal components to the direction of dominant compliance. At positive angles, the projections of the force response from the tangential and normal components act in the same direction and at negative angles in the opposite direction. In the case of tilt angle as a limiting variable in SLD in Fig. 8b, a wider region of instability for negative values is noticeable. This can be explained by the smaller mean distance of the engagement points from the axis of rotation (the cutting edge spends more time in the engagement) for negative angles, which leads to a larger force response of the process and therefore lower stability. In contrast to the lobe diagram on a cylindrical cutter, where the number of teeth only scales the limiting axial depth of cut, the non-linearity of the stability problem for a ball cutter is clearly evident from the different shape of the lobes for different numbers of tool teeth, see Fig. 9.

Discussion
The proposed SLD calculation method builds on Otto's approach and extends it to a cutting process with nonlinear behaviour of the cutting force gradient on selected limiting parameters, such as ball milling. It also generalizes the approach of Ozturk and Budak [7] and Ozturk et al. [8] to calculate the Jacobian of the cutting force (direction matrix) to a general tool shape and formulates it in an analytical form instead of a discretized one. The analytical formulation has a major advantage besides speed of calculation in the ability to calculate the derivatives of the Jacobian (directional matrix) with respect to the process parameters.
This allows, first of all, to efficiently find the solution of the limiting parameter in the SLD, but also to determine which of the process parameters has the dominant effect on the stability of the process at a given point in the SLD. This approach is demonstrated by calculating the SLD on a ballend milling where, compared to previous articles, the angles defining the engagement have been refined by analytically describing the surfaces defining the tool envelope,  the unmachined surface of the workpiece and the surface formed by the previous groove. The angles of engagement were determined from the equations defining the intersections of these surfaces. This is a more reliable method in 3D geometry than approaches based on difficult to prove reductions to 2D geometry and trigonometry. The engagement areas required to calculate the Jacobian for ball milling are limited to the case where only the ball portion of the tool is in the cut.
The potential of this approach lies in its use in simulation software for general tools, where the calculation of the surface integral can be performed by numerical quadrature while taking advantage of the existing discretization of the removed material.
The cutting force Jacobian depends nonlinearly on all relevant technological parameters which led to a need to work out a more general condition on machining stability limits than derived by Otto et al. [21]. This condition is applicable not only to ball-end milling but also to some other end milling tools (bull-end, toroidal) or to milling tools with round inserts or viper corner geometries, where the geometry leads to a non-linear Jacobian with respect to the depth of cut. It can also be applied to face or peripheral milling, whose Jacobian (directional matrices) is linear in the axial depth of cut but non-linear in the radial depth of cut. The calculation of the limits in case of the non-linear problem is easier when the derivatives of the Jacobian with respect to the sought limit parameter are known. This is the case for the presented case of ball-end milling.
The method of calculating the mean value over tool revolution through the surface integral is used in the paper to calculate the mean cutting force Jacobian, but this identity can be straightforwardly modified for the calculation of the mean cutting force itself and used, for example, to identify cutting coefficients in generally oriented tool in a similar way as Gradišek et al. [31] did for a end mills with no inclination.

Conclusion
An alternative mathematical formulation of the ZOA-based machining stability prediction was developed for ball-end milling. The method is suitable for complex engagement geometries and operations where calculation of the directional matrix is difficult and the previous solutions relied on material removal software for tool engagement determination. Another advantage of the method of calculating the Jacobian (directional matrix) using the surface integral compared to its calculation from the definition of the ZOA method is that it is not necessary to know the time evolution of the Jacobian and therefore significantly fewer integration points are sufficient, especially if the area of the engagement is small compared to the area of the tool package.
The proposed calculation of the mean cutting force Jacobian was validated independently in material removal simulation software. The results of the proposed method were compared with averaged time-dependent cutting force Jacobian, which was obtained numerically from cutting forces calculated on toolpaths sequentially perturbed in the feed, cross-feed and surface normal directions.
The cutting force Jacobian matrix (directional matrix) obtained by the proposed method is used in a generalized stability condition for the delayed system extending the approach developed by Otto for non-linear eigenvalue problems. The analyticity of the formulation enables differentiation with respect to process parameters such as tilt, lead angles, or radial and axial depth of cut, which could be useful for machining stability optimization.
where s min , s max denote lower and upper limit of cutting edge parametrization. The merit of the formulation is that while the previous procedure works with a parametrization based on tool rotation and a curve parameter for the cutting edge, the reformulation into a surface integral no longer depends on the surface parametrization and it is therefore possible to choose the parametrization in which the calculation is the simplest. The derivation of the formulation is based on a textbook formula for a surface element in the curvilinear coordinate system, which is The cross product can be conveniently written as which leads to formulation used in the integral above where ̂ is the tool axis direction unit vector.

A.4 Script for cutting force Jacobian calculation
The calculation of the Jacobian for milling spherical ends leads to matrices that are too complex to be reasonably presented in text form. The code below in Wolfram cos sin 2 (tan L cos − tan T sin ) + tan L cos 2 + sin sin cos tan T � sin 2 sin 2 + cos 2 � + cos sin (cos − tan L sin sin ) sin cos (tan T cos + tan L sin ) + sin 2 ⎞ ⎟ ⎟ ⎟ ⎠ . Mathematica follows the calculation described in the paper. The division into the functions shown is made with readability in relation to the article in mind, not code efficiency. For faster computation, it is preferable to analytically precalculate some of the functions (e.g. matrix [T]).