Chatter-free milling strategy of a slender Blisk blade via stock distribution optimization and continuous spindle speed change

The machining process of Blisk blades poses multiple challenges due to high requirements on surface quality and precision combined with high dynamic compliance of the thin-walled blades. Avoidance of chatter is thus of high priority in Blisk blade machining. However, the geometry of the Blisk blade array where the tool must fit between individual blades significantly limits the possibilities of controlling stability through the relative orientation of the tool and workpiece. Thus, the main parameters that can be used to control the stability of the process are the distribution of stock allowance and the spindle speed. Due to the effect of material removal on the blade’s dynamic properties, spindle speed needs to be adjusted throughout the machining process to keep it within the continuously changing stability gaps. This paper describes in detail an optimization procedure for the design of stock allowance distribution in such a way that a continuous spindle speed modulation is possible that avoids chatter throughout the machining process by maintaining spindle speeds within stability gaps. The presented algorithm uses finite element analysis software to simulate the effects of stock allowance distribution and material removal on workpiece dynamical properties. This information is then coupled with a stability model based on the Jacobian of the cutting force with respect to the regenerative deflection to identify the varying stability gaps throughout the machining process. The proposed method was experimentally verified.


Introduction
Flexible thin-walled workpieces such as Blisk blades are significantly represented in key sectors of the aerospace, automotive, and energy industries. The low static and dynamic stiffness of these products leads to an increased risk of deformation and instability (vibration) caused by the force interaction between the tool and the workpiece. At the same time, the surface and precision requirements on the end products are high. These physical constraints and production requirements make efficient machining of these parts a challenging task, and, in many cases, manual finishing of such workpieces is necessary. This significantly increases the production time and the risk of additional geometric inaccuracies. The most serious dynamical behavior issue is the selfexcited vibration called chatter, which is caused by the force interaction between the displacement of a compliant structure and the force response of the cutting process to this displacement and its evolution, which is recorded on the surface of the workpiece. The positive feedback that leads to oscillation is due to the fact that the stiffness of the process depends on the regeneration of chip thickness, an effect first considered by Tlustý and Špaček [1] for turning. In the case of milling, the problem is more difficult because of the time-dependent force response of the cutting process. To simplify this problem, Budak and Altintas [2] developed an approximate method to solve the vibration stability for cylindrical end mills, which is called zero-order approximation (ZOA). The method is based on averaging the timedependent directional cutting force matrix, i.e., it uses its zero-order Fourier series term. It was compared with more accurate time-domain simulations and found to be reasonably accurate for tools with standard geometry (tools without serrations or variable helix angle) and operation (constant spindle speed) [3]. For a recent state-of-the-art overview on chatter and geometric error in thin-walled machining, see [4,5].
Methods for calculating milling stability were first developed for cylindrical tools in the context of 3-axis milling due to the simpler geometry of the process, and some authors generalize those approaches to 5-axis milling. However, these methods can only be used for some specific cases of 5-axis milling [6][7][8][9]. There are two major differences between these specific 5-axis cylindrical milling operations and the general five-axis ball milling that complicate stability calculations. First, for cylindrical tools, the directional of the cutting force is directly proportional to the depth of cut (the directional matrix does not change with depth of cut). No such simple relationship exists in the case of general 5-axis ball-end milling. The second difference lies in the detection of cutter-workpiece engagement region. In the specific cases of cylindrical milling, the engagement region can be evaluated using simple trigonometric formulas with the use of a 2D representation of the cutting process. In contrast, no such direct method of tool engagement region exists for general ball-end 5-axis milling. For this reason, it is often advantageous to use simulation software (such as in [10] or [11]) to calculate tool and workpiece engagement.
Three main approaches to computation of cutter-workpiece engagement regions can be identified in the literature: analytical, solid, and discrete modeling methods.
The main advantage of analytical modeling methods is their precision and low computational requirements (compared to the solid and discrete methods). However, the concrete analytical methods are applicable only to specific settings (three/five-axis machining, point-end/flank milling, specific types of cutters), thus limiting their application scope. In [12] an analytical method of calculating the cutterworkpiece engagement region for toroidal and flat-end cutter during semi-finish milling of sculptured parts is presented. In [13], the authors present an analytical model of CWE and instantaneous chip thickness in three and five-axis milling with flat workpiece surfaces. An analytical CWE computation method based on the integration of the arc-surface is proposed in [14]. In [15], the author presents an analytical model of the cross-sectional profile of a groove produced by helical milling with a ball-end mill.
Solid modeling methods rely on the exact representation of the tool and workpiece geometry by a hierarchy of interconnected geometrical primitives called the boundary representation (B-rep). These methods calculate the CWE by applying Boolean operations to the workpiece and tool swept volumes and surfaces. Solid modeling methods are generally very accurate, but suffer from high computational requirements. Furthermore, these methods often require the utilization of commercial geometric modeling kernels. For some of the applications of solid modeling methods to CWE computation, see, e.g., [10,[16][17][18][19][20].
Discrete modeling methods aim to simplify complex solid modeling calculations while retaining required precision by applying surface and volume representation methods that typically originate in the computer graphics domain. Most commonly applied variants of these methods include the Z-buffer [21][22][23] and distance field representations (see, e.g., [24][25][26][27][28]).
The method used to calculate CWE in this paper relies on a distance field representation of workpiece material and is essentially the same as that in [24] and [25]. There, the CWE was computed by defining a sampling grid on the tool's surface and then identifying the grid points that are in contact with the workpiece material (while taking into account the instantaneous feed direction to correctly identify the egress points). For a description of the CWE calculation applied in this paper, see Appendix C).
In some cases, high surface roughness and low accuracy caused by chatter can be avoided by using special tools. The effectiveness of chatter suppression of standard, variablepitch, and crest-cut tools in milling thin-walled workpieces is compared in the recent article [29]. Although the crest-cut geometry compares favorably with the other types of geometry in this regard, its application is limited in point milling, where the engagement surface is small compared to flank milling. Another possibility is to change the lead and tilt angle settings [30,31]. In some cases, however (such as with monolithic Blisk wheels), the combination of tool and workpiece geometries induces severe restrictions on the permissible lead/tilt angle ranges. Therefore, alternative strategies are required in order to eliminate vibrations and decrease surface roughness. These strategies include spindle speed control [32] and stock allowance optimization coupled with spindle speed control [33][34][35]. The latter articles took into account the change in structural dynamics of thin-walled blades relative to tool position when optimizing the feed and spindle speed distribution to reduce the risk of chatter. To do this, a structural modification algorithm based on the FE mesh of the final workpiece was applied in [34]. In [35], the shape of the stock allowance distribution along the workpiece surface is modeled by a combination of harmonic and exponential functions and, together with tool-workpiece orientation, optimized with respect to machining stability.
In structural dynamics modeling of thin-walled workpieces, it is not only the low compliance of the workpiece that complicates machining but also the fact that the workpiece compliance changes significantly depending on the position of the tool [36,37] and with respect to the distribution and continuous removal of material from the workpiece. When considering the effect of stock allowance distribution on the workpiece surface, it is necessary to balance two opposing effects -on the one hand, higher allowance stiffens the thin-walled workpiece but on the other hand, it increases the force response of the cutting process. The variable dynamical properties limit the use of the otherwise well-established modal hammer tests for this type of task, where FEA analyses described in [38,39] are used to determine the dynamic properties of the workpiece at the cutting location. In [40] the stability of a machining system is modeled during high-speed milling, taking into account the dynamic properties with respect to both the tool and the actual position on the workpiece. The stability prediction in 5-axis milling of a thin-walled blade with a predetermined stock allowance based on the FE model and the geometric simulation of the process in the time domain is also addressed in [41]. The authors of [42] focused on the effect of self-excited vibrations on the development of surface roughness in thin-wall milling. In [43] the variation of the FRF with respect to material removal was approximated using time-efficient reduction methods, combining dynamic substructuring and matrix perturbation methods. The finite element method in combination with the structural dynamic modification method was used in [44] to determine the changes in dynamics of a curved surface workpiece. The authors of [45] used FEA with both shell and cubic elements to capture the time-dependent dynamics and developed the resulting mass and stiffness matrices into a Taylor series. Recently, the authors of [46] developed a hybrid method of simulation of position dependent dynamics using FEA with shell elements. This method uses the shell element formulation to update workpiece dynamics by applying the matrix perturbation method. The accuracy of the results is further improved by using either full FEA solutions or experimentally obtained FRFs to reduce the accumulated perturbation errors.
In the proposed method, a parametrized FEA model was applied to iteratively evaluate the dynamics of the workpiece throughout the machining process. While computationally more demanding when compared to the above-mentioned methods of matrix perturbation and Taylor series expansion, it proved to be easier to implement and its evaluation was not prohibitively time-consuming. The time requirements were further lowered by simplification of the material removal process and by avoiding some remeshing procedures.
This article presents a chatter avoidance milling strategy of a Blisk blade that considers both the semi-finish and the finish phases. Specifically, the presented optimization strategy focuses on finding a distribution of stock allowance of the semi-finish end product that provides increased workpiece stiffness, while simultaneously allowing for a continuous spindle speed change during the finish phase that suppresses chatter. To achieve this, the effect of stock allowance distribution and its removal on dynamical properties of the workpiece is simulated via a finite element (FE) model with parametrically controlled stock allowance distribution. The dynamical properties are then combined with the force response of the cutting process to identify feasible stability gaps. Similarly to the previously mentioned stock allowance optimization methods, the proposed approach balances the effect of increased workpiece stiffness due to added stock allowance against the increase in force response of the cutting process associated with a higher depth of cut. A backtracking algorithm is applied to find the distribution of stock allowance for which the spindle speed can be continuously modulated throughout the machining process so that its value always lies within the (continuously changing) stability gap. These computed spindle speed values are then used to modify the process NC code.
The main novelty of the proposed method thus lies in maintaining continuous spindle speed throughout the finishing process in such a way that spindle speed values never leave the highly variable stability gaps. This is achieved by considering the effect of continuous material removal on workpiece dynamics and iteratively searching the solution space for those stock allowance distributions for which such continuous spindle speed modulation is possible.
This article is organized as follows. Section 2 presents an outline of the optimization process presented. Section 3 presents the workpiece description and discusses the chosen milling strategy. Section 4 discusses the general design of the applied material removal simulation. Section 5 describes the structure and parameters of the finite element model used to simulate the removal of stock allowance removal and to evaluate the dynamical properties of the workpiece throughout the optimization process. Section 6 describes the combined stock allowance and stable spindle speed optimization process along with a description of the applied cutting force model. Section 7 focuses on the methods of calculation of the speed function and subsequent modification of the NC code. Section 8 presents the experimental verification. Finally, Sect. 10 summarizes the benefits of the presented optimization approach and discusses future research possibilities. The article is concluded with three appendices, Appendix A containing details on the resulting stock allowance distribution obtained in the course of experimental verification of the presented method, Appendix B containing detailed photographs of final workpieces produced by the machining strategies compared in this study, and Appendix C describing the material removal simulation software used for tool engagement region identification.

Outline of the optimization strategy
This section presents an outline of the proposed optimization process for the reader's convenience. A more detailed description of the individual optimization steps is provided in subsequent sections.
In the finishing phase of the machining process, flexible Blisk blades are typically machined from top to root with a helix toolpath strategy separated into several levels. A different set of constant stable cutting conditions is specified at each level based on experimental results and previous experience. However, this approach leaves undesirable traces on the workpiece surface at level transitions and can lead to suboptimal machining times. Moreover, the dynamics of Blisk blades continuously change during the machining process. This, coupled with the high compliance of the blade, significantly increases the risk of chatter.
In an effort to mitigate these disadvantages, this paper presents an extension of this standard machining process which includes optimization of the stock allowance distribution as well as continuous spindle speed control resulting in a single uninterrupted chatter-free machining process (see Fig. 1). In order to evaluate the influence of distribution of stock allowance on the stability of the finishing process, a method of material removal simulation is required. In the proposed method, the gradual removal of stock allowance was simulated by discrete removal of stock allowance in several mutually parallel horizontal sections (denoted as strips in the following), see discussion in Sect. 4.
The modeling of the material removal itself and its effect on the blade dynamical properties was performed with the use of a finite element (FE) model designed in the ANSYS software. In this model, the finished workpiece volume was modeled by quadratic solid elements, while the stock allowance was modeled by shell elements of adjustable thickness (see Sect. 5.1).
In order to evaluate the effects of stock allowance distribution on the in-process changes in workpiece dynamical properties, several points were defined on each considered strip of material. The effects of material distribution on the stability of the machining process within any given strip were evaluated by combining the local stability properties obtained at these points (see Sects. 5.2 and 6.4).
In this way, the stability limits a max p i (thickness of stock allowance) and the corresponding stable spindle speed intervals I max i were computed (see Sect. 6) for any simulated stock allowance distribution and any specific point during the machining process. Since the dynamical properties of the workpiece at any point during the machining process are influenced by the material distribution on the not yet machined parts of the workpiece, the final stock allowance distribution was computed in reverse order to the actual machining process starting at the last strip to be machined and ending at the first (topmost) strip.
The algorithm responsible for the computation of the final stock allowance distribution was a backtracking algorithm iteratively searching the space of possible stock allowance distributions for those stock allowance distributions for which stable spindle speeds could be continuously modulated during the whole machining process.
After such a stock allowance distribution (and its associated combination of spindle speed intervals) was found, a spindle speed function was computed that represented the desired continuous change of spindle speed throughout the whole machining process (Sect. 7.1). The spindle speed function was then used to modify the NC code by extending the NC code commands with the information about the respective spindle speed values (see Sect. 7.2).
For a schematic representation of the optimization process described in this section, see

Blisk blade -characteristics and milling strategy
A segment of one Blisk bucket from a monolithic Blisk wheel was considered for the purposes of the experiment. The workpiece material was an aluminum alloy EN AW 7075. The dimensions and shape of the blade are shown in Fig. 3. For the mechanical properties of the alloy used, see Sect. 5.1, Table 1. The helix machining strategy with stepover 0.4 mm was applied (see Fig. 4). The input parameters for the optimization method were the range of depth of cut 0.1-0.9 mm and the range of spindle speed 10000-15000 rpm . The range of depth of cut was specified on the basis of previous experimental experience, taking into account material properties, tool geometry, and cutting conditions. The limits of the spindle speed range were based on the specifications of the spindle used in the experiment. The admissible range for the lead angle was 80-85 • and the range for the tilt angle was 0-4 • . The lead and tilt angle ranges were forced by the geometry of the Blisk wheel (to avoid tool collision with neighboring blades during the machining process). To further reduce the risk of chatter, the tilt angle closest to the normal to the machined surface (i.e., 80 • ) was chosen. The value of the lead angle was set to 4 • in order to avoid zero cutting speed when machining the radius of the blade root where the tool tip can come into contact with material.

Method of material removal simulation
Taking into account the chosen milling strategy and in order to reduce the overall complexity of the optimization process, the surplus material was modeled as a series of mutually parallel horizontal sections (strips) covering the whole surface of the Blisk blade (i.e., both its suction and pressure sides). Furthermore, a regular grid of points in which local dynamic properties were evaluated was defined on the pressure side of the Blisk blade (see Fig. 5).
The continuous process of removal of stock allowance was thus simulated (discretely) by removal of individual material strips, neglecting the effect of gradual removal of material in the horizontal direction. This decision was motivated by the fact that for finishing and semi-finishing operations with spot machining on the selected knife, the effect of material removal within a single line is negligible. The mass of material removed per row is on the order of lower tenths of a gram, and the effect on stiffness is also marginal given the chosen strategy of gradual removal of material from the top of the blade to the root.

Finite element model of the Blisk blade
This section presents details of the construction of the finite element model of the blade (Sect. 5.1), the methods of evaluation of workpiece dynamical properties (Sect. 5.2) and model verification methods and results (Sect. 5.3).

FE model of slender blade
The finite element model of the blade was created in ANSYS 2021 R2 using a combination of solid and shell elements. The final blade geometry (including both leading and trailing edge radii) was represented by a quadratic mesh of solid elements. The main thin body of the blade (excluding the root radius and base) was swept with cubic elements (SOLID186) with 48 divisions (division size 1.8 mm and element size 1.35 mm ). Blade's base and root radius were modeled using a free tetrahedron mesh (SOLID187) with element sizes ranging from 3 to 8 mm.  The stock allowance was modeled using a layer of quadratic shell elements (SHELL281) coated on both the suction and pressure side (including root radius), creating a conformal mesh with solid elements (see Fig. 6).
To model the distribution of stock allowance, the blade surface was divided into 12 individual strips according to a division of the solid swept mesh such that the width of each strip spanned four divisions of elements. The stock allowance distribution was controlled parametrically by setting the shell thickness (measured in the normal direction to the final geometry) for each of the 12 strips (see Figs. 5 and 6). Strips with desired zero stock allowance thickness (i.e., already machined strips) were modeled by setting their thickness to an auxiliary value of 1 ⋅ 10 −5 mm . This value was small enough to not influence the precision of the model while preserving the element mesh of the model (which would have to be recalculated if the shell elements were disabled). The degrees of freedom at the circular surface of the model base were removed.
The mechanical properties of the blade material are given in Table 1.

Dynamical compliance modeling and evaluation
The blade's dynamical properties were evaluated on a regular 12 × 9 grid of points. This grid consisted of selected nodes of the swept grid on the pressure side of the blade. Each of the 12 rows of 9 grid points lied in the middle of the corresponding strip (see Fig. 5). In the following, the number of grid columns and rows shall be denoted as N col and N row , respectively. A local coordinate system (LCS) was defined in each of the grid points P i,j as a matrix of orthonormal column vectors where the tangent vectors t i,j were defined as and i.e., as normalized direction vectors connecting two consecutive grid points within a strip. This approximation was based on the planned toolpath design (see Fig. 4) in which the tool moves in a horizontal direction when machining the blade sides. The normal vectors n i,j were set as surface normal unit vectors obtained from an STL model of the finished workpiece exported from the ANSYS software. Finally, the binormal vectors were defined so as to convert LCS P i,j into right-oriented orthonormal systems.
The differences in dynamical behavior between the suction and pressure sides were ignored (hence the location of the dynamical grid points on the pressure side only). This simplification was justified by the slenderness of the blade

Model verification
This section describes the results of two separate experiments that aimed to verify the accuracy of the material addition/removal simulation method (Sect. 5.3.1) and the accuracy of the FE model of the finished blade geometry (Sect. 5.3.2) from the point of view of eigenanalysis.

Comparison of solid versus solid+shell models
The method of simulation of material addition using the solid+shell combination was validated by a comparison of two dynamical models, both representing the final blade geometry with stock allowance distribution of constant thickness of 0.5 mm . In the first model, the target geometry was completely meshed using solid elements only. In the second model, the combination of solid and shell elements described in Sect. 5.1 was applied. To compare the two models, the frequency response functions (FRFs) were evaluated at the most compliant point at the top of the blade. The comparison (see Fig. 7) revealed minimal differences between the two models. Specifically, the largest relative frequency difference in the first 20 eigenmodes was evaluated at 0.32% corresponding to an absolute difference of 41.69 Hz at the eight eigenfrequency 13082.98 Hz of the solid+shell model.

Comparison of the finished geometry FE model versus the finished workpiece
To verify the accuracy of the FE model of the final geometry (i.e., with shell thicknesses set to the value of 1 ⋅ 10 −5 mm ), the reduced-order state-space model with a constant damping ratio of 0.35% was compared with the modal measurements of the finished workpiece. As in Sect. 5.3.1, the FRFs at the most compliant point of the blade were compared (see Fig. 8).
The modal measurement was realized by clamping the blade in a collet and performing a tap test using an impact hammer of very small mass ( 4.8 g ) to avoid double taps. The response was measured with a non-contact Laser Doppler Vibrometer to avoid the mass impact of traditionally used attached transducers. The experimental equipment and settings are summarized in Table 2.
The modal measurement was performed to verify that the FE model provides a sufficiently precise frequency agreement when compared to the finished workpiece. The results show a very good (less than 5% relative difference) frequency agreement at the first four eigenmodes (see Table 3 and Fig. 8). Moreover, the frequency difference was less than 1% at the first eigenmode.
In addition, the damping ratio at each eigenmode was evaluated from the tap test results as an average of two methods: the Least Squares Complex Exponential (LSCE) and the Ibrahim Time Domain (ITD). Both these methods operate in the time domain and are therefore suitable for structures with low damping. The first eigenmode exhibited  an especially low damping of 0.04% that led to a sharp peak of dynamic compliance with maximal value of 5.6 N∕mm (see Fig. 8). The relatively large amplitude differences in frequency responses at the eigenmodes were not problematic, since the optimization method presented primarily focuses on avoiding the eigenmode frequencies which were matched by the FE model with sufficient accuracy.

Stock allowance distribution and spindle speed optimization
This section describes in detail the individual components of the optimization process that lead to the final distribution of stock allowance and the associated stable speed intervals. The discussed components are: the details of stock allowance thickness and spindle speed parameter discretization (Sect. 6.1), the stock allowance distribution optimization (Sect. 6.2), the stability evaluation methods (Sect. 6.3) and the discrete lobe diagram computation (Sect. 6.4).

Stock allowance thickness and spindle speed discretization
Taking into account both the limitation of the technological process that does not allow for machining with arbitrary precision and the computational requirements of the optimization process, the interval of possible material thicknesses and the interval of achievable spindle speeds were discretized. The vector of the material thicknesses considered was selected as and the interval of the spindle speeds considered was selected as It should be noted, however, that for the purposes of the evaluation of the stable spindle speed interval on the 12-th (topmost) strip, a finer sampling of a narrower spindle speed interval was used.
The finer sampling of the frequency value vector S p ⋆ on the 12-th strip was forced by the higher dynamical compliance of the blade near the top of the blade.

Optimization of stock allowance distribution
The main aim of the stock allowance distribution optimization algorithm was to define a sequence of stability limit values {ap max i } 12 i=1 and associated stable spindle speed intervals for all material strips under the constraints and The constraint (1), although not strictly necessary, was primarily introduced in order to limit the size of the search space of the considered stock allowance distributions.
The constraint (2) was introduced in order to guarantee that a continuous transition of stable spindle speed was possible across boundaries of adjacent material strips. The parameter Δ was set to 50 rpm for all strips, except for the boundary between strips 11 and 12 where its value was set to 25 rpm . This change of value was forced by the higher dynamical compliance of the blade near the top.
When computing the stability limit ap max i on the i-th strip, the FE model was evaluated for all material distributions of the form The H i,k parameter combinations correspond to all stock allowance distributions at the point of machining of the i-th strip that the model considers (i.e., at the point of the finish phase when the material has already been removed from all the strips with higher index). In order to efficiently evaluate these material distributions, it is beneficial to parallelize their computation with respect to the discrete material thickness values in A p . This was achieved by executing the  ANSYS script from MATLAB using the system command inside a parallel for loop (parfor command).
In the proposed algorithm, the stable spindle speed limit was selected as the maximal thickness of the stock allowance compliant with the constraint (1) for which there existed a stable spindle speed interval respecting the constraint (2). If there was more than one stability interval that satisfied (2), the stability interval of the greatest length was chosen.
Of course, any particular choice of ap max i and I max i on the i-th strip has an influence on the size of the solution space in subsequent (higher) strips. For some choices, no solution can be found at subsequent strips. For this reason, the algorithm included a backtracking feature. If at any strip a solution compliant with constraints (1) and (2) could not be found, the solution at the previous strip was modified by first considering alternative choices (if any) of stability intervals or else considering a lower value of ap max i at the previous strips. The described optimization method presents a considerable degree of freedom when selecting the stability limit and stable spindle speed interval on the first strip (i.e., at the root of the blade). A reasonable approach is to select the highest stability limit for which there exists a stable speed interval of some minimal length. In the specific case described in this paper, the process was deemed stable by the optimization method for all combinations considered of A p (stock allowance thickness values) and S p (spindle speed values) on the first strip. Therefore, the selected values on the first strip were and

Stability evaluation
This section briefly describes the applied cutting force model along with the method of cutter-workpiece engagement calculation (Sect. 6.3.1)and stability determination (Sect. 6.3.2).
The stability of the machining process is determined by the coupling between the force action of the structure and the process. In the case studied, the workpiece compliance was dominant compared to that of the tool and was determined by the FEM of the workpiece, which had damping tuned according to the experimentally determined FRFs.

Calculation of cutting force and its Jacobian
The modeling of the process force effects is based on the approach of modeling the force between the tool and the workpiece by integrating the specific force on the cutting edge element along a curve representing the cutting edge [47]. The force acting on an element of cutting edge of width dw can be expressed as where g is a characteristic function whose value is 1 if a point s on the j-th cutting edge rotated by the angle is engaged and 0, if it is not. The vectors K e , K c denote the experimentally identified edge and cutting force coefficients, respectively, and h is the dynamic undeformed chip thickness approximated by the Martellotti formula [48] where f t denotes feed per tooth vector, denotes the tool and workpiece regenerative dynamic displacement calculated from tool-workpiece dynamic displacement , and denotes time delay between two consecutive cuts. The elemental chip width dw is defined as the length of the cutting edge element projected onto a plane perpendicular to the tangential direction. For an example of a local coordinate system basis at two different points, see Fig. 9.
The specific cutting force coefficients were determined experimentally for the EN AW 7075 aluminum alloy and a helix angle of 30 • , which corresponds to the helix angle on the shank of the ball-end tool used for the machining. The values of the cutting coefficients were identified using a method by Gradišek et al. [49]. In this method, the cutting coefficients are estimated from mean cutting forces computed separately for multiple different values of feed per tooth (in this case, 0.08-0.2 mm ). These mean cutting force values were computed using tool engagement information obtained from the material removal software MillVis (see Appendix C. The setup of the cutting force coefficient identification experiment is presented in Fig. 10.
The identified coefficient values were: where the vector elements correspond to coefficients in the tangential, normal and binormal directions, respectively.
In the presented analysis, we consider a steady (i.e., quasi-stationary) state, where the equations describing the forced oscillation and the self-excited oscillation can be separated on the basis of Floquet theory. The Jacobian of the cutting force with respect to the regenerative deflection is a crucial component of the machining stability formulation (see [50]).
The cutting force Jacobian averaged over one revolution necessary for the zeroth-order approximation based stability lobe diagram calculation can be expressed as In order to evaluate the above expression, it is necessary to be able to determine whether the blade is in engagement for relevant values of s, i, and . Formally, this information is provided by the characteristic function g. A description of its evaluation is briefly discussed in the following section.
As mentioned in the Introduction, it is difficult to determine the engagement region for a general tool and toolworkpiece orientation. A common method (see [51,52]) for determining the engagement conditions is to use a material removal simulation software that generates a matrix of zeros and ones determining whether grid points on the tool surface are in engagement for a discretized tool rotation and blade position. This method was applied in the presented study. The material removal simulation used is shown in Fig. 11. A brief description of the software can be found in Appendix C.

Stability exponent calculation
Averaging the cutting force gradient leads to a system of delayed differential equations with constant coefficients. The dynamic behavior of the system under a dynamic load acting on the workpiece is mathematically described in the Laplace domain as where is a mass-normalized matrix of eigenvectors and U is a generalized eigenvector matrix. The sign of the real part of the parameter s determines if the machining process is stable or unstable. This problem was solved using the iterative method described in [53] (see also [54] for a more detailed description).

Discrete lobe diagram computation
After dynamical properties of the i-th strip grid points are computed for all material distributions of the form described in (3), a discrete version of the classical lobe diagram is computed for every grid point P i,j of the i-th strip. This discrete lobe diagram is defined as a N ap × N sp matrix L i j = {L i j (k, l)} (where N ap and N sp denote the lengths of the A p and S p vectors, respectively), with individual elements defined by The condition C i j (k, l) in (4) is satisfied whenever (5) holds: The evaluation of condition C i j (k, l) depends on several variables (lead and tilt angles, tool geometry, local surface geometry). The evaluation method is described in detail in Sect. 6.3.
After the discrete lobe diagrams of the individual points are evaluated, a Discrete Lobe Diagram Union (DLDU) matrix L i = {L i (m, n)} for the i-th strip is defined as It is immediately seen from (4)-(6) that an element of L i is zero precisely when the machining process is stable for l) is satisfied, 1 otherwise.

(5)
The cutting process is stable at point P i,j when removing surplus material of thickness Ap(k) from the i-th strip at spindle speed S p (l) given surplus material distribution H i,k .
the specific combination of conditions for every grid point of the i-th row. Therefore, the intervals of spindle speeds for which machining is stable along the entire strip can be identified from L i (while discounting the effects of gradual material removal in the horizontal direction) for every stock allowance thickness value in vector A p by examining the set of zero values of the respective row of L i . A visualization of the composition of the DLDU from discrete lobe diagrams of individual points is presented below (Fig. 12).

Application of optimization results to the machining process
This section discusses how the results obtained by the methods described in Sect. 6 can be applied to the machining process, i.e., how to associate individual NC code position commands with their respective spindle speed values in order to achieve a chatter-free finish process.

Spindle speed function
After the final distribution of stock allowance thicknesses along with the respective stable spindle speed intervals (height of the rectangle) and stability interval I 6 (width of the rectangle) is found, it is necessary to construct a function Ω , that defines the dependence of spindle speed on the tool contact point coordinates. Specifically, the function is defined to depend on the Z-coordinate of the tool contact point, as the chosen machining strategy consists of gradual removal of material in the Z-axis direction.
The function Ω is expected to simultaneously satisfy the following requirements: ) is necessary in order to simplify the machining process and avoid the issue of surface deterioration at the points of spindle speed discontinuity. The requirements on higher order continuity (R2.) and minimal oscillatory behavior (R3.), while not strictly necessary, are introduced to ensure that the spindle speed variation prescribed by Ω is within the kinematic capabilities of the spindle speed control.
In order to satisfy the requirements [R1.]-[R3.], the function Ω is computed as a result of the optimization problem (8). This optimization problem aims to minimize the total variation of Ω . The function Ω could be defined as piecewise linear, but in order to guarantee a higher degree of continuity, Ω was defined as a shape-preserving piecewise-cubic Hermite interpolation spline (PCHIP) instead. PCHIP functions display a higher degree of continuity while simultaneously preserving shape and monotonicity of interpolated data (for a more detailed description of monotone spline interpolation, see [55,56]). This makes them a better alternative to either piecewise linear interpolants or standard spline interpolants (which typically display a higher degree of oscillation).
The function Ω is thus defined as a shape-preserving piecewise-cubic Hermite interpolation spline that interpolates points where the coordinates z k are defined as: The values of z k defined in (7) correspond to z-axis coordinates of the bottom, middle, and top of the individual material strips, respectively.
The function values s k are obtained as the results of the non-linear optimization problem The boundary conditions specify that the values of Ω at the bottom of the lowest strip and in the middle of each strip must lie within the respective intervals of stable spindle speeds (Eqs. (9) and (10), respectively). The values of Ω at the common boundaries of two consecutive strips must lie at the intersection of the respective stable spindle speed intervals (11) and the value at the top of the top strip is fixed in the middle of the corresponding stable spindle speed interval (12).
The spindle speed function used in the experiment presented in this paper is displayed in Fig. 13.

NC code modification
CNC milling centers are typically programmed with a safety feature that defines a spindle speed tolerance interval centered around the (constant) programmed spindle speed. If the spindle speed exceeds this interval at any time during the machining process, the movement of the machine axes is stopped and a warning is raised. This safety feature is too strict to allow for a continuous spindle speed change during the machining process. Therefore, this limitation was bypassed on the MCU700 machine used in experimental verification by a custom modification of its PLC programming. This modification extended the safety feature, allowing the machine tool to receive new programmed spindle speed commands during the machining process without interrupting axes movement.
In order to produce the final NC code, the tool contact point coordinates were first calculated for each line of the unmodified NC code. These coordinates (or rather their Z-axis values) were then used to evaluate the optimized

Experimental validation
The proposed machining strategy was verified by a series of experiments on a MCU700 5AX five-axis milling machine tool (see Fig. 14) with a Heidenhain TNC640 machine control system. In all experiments, the ISCAR MM EB100A07-4T06 ball-end mill with 5-mm radius was used. A total of four Blisk blades, denoted A, B, C, and D, were machined in the course of the experiment. Blades A-C were machined using conventional machining strategies with (mutually different) constant spindle speeds and constant stock allowance thicknesses, while the D blade was machined using the optimized strategy described in this paper. The strategies of A-C were the result of several (unsuccessful) attempts to eliminate chatter through a variation of the cutting conditions based on the extensive practical experience of the machine operator. For a comparison of the parameter settings used in machining the test blades, see Table 4. In order to compare the surface quality of the experimental results, four measurement areas were selected on the blades surfaces, two on each side. The dimensions and  positions of these areas are presented in Fig. 15. The roughness parameter Ra was calculated on each area based on five separate measurements taken in the direction of the toolpath (denoted by the purple arrows in Fig. 15). The results of the roughness measurements are presented in Fig. 16. These results show that of all the tested configurations, only the D blade machined by the proposed strategy satisfied this threshold in all measured areas.
The locations of the individual areas were selected based on the following reasons. Area 1 was selected because it covers the part of the surface on the suction side with the highest tendency for surface deterioration due to machining instability. Area 2 corresponds to a rigid part of the surface where no chatter was observed for any of the strategies tested. The surface roughness value measured on the area 2 therefore served to verify that the prescribed roughness limit of 0.8 μm was achievable with the chosen machining strategies (compare with the results for area 2 in Fig. 16). The location of Area 3 was chosen because it covers the part of the surface on the pressure side with the highest tendency for surface deterioration due to machining instability. Finally, area 4 covers the transition area from the highly flexible part near the top to the more rigid structure near the root of the blade.
The roughness values in areas 1 and 3 confirm the results of the visual inspection (see Fig. 17 and Appendix B for the comparison of blade surfaces) where significant undercuts and deterioration in surface quality can be observed in these areas for the configurations A-C. Specifically, compare the surface details of area 1 presented in Fig. 18. The roughness values obtained at the most rigid measurement area, area 2, confirm that (excluding the effects of machining instability) the cutting conditions of the machining strategies were chosen correctly with respect to the specified roughness threshold. The comparison of surface roughness values measured in area 4 shows that the effects of machining instability were not prominent in this area (all strategies except strategy B fall below the roughness threshold of 0.8 μm ). This can be explained by the increased rigidity of the blade in this area. Fig. 18 Comparison of surface detail of the tested machining strategies. Images refer to the measurement area 1 (see Fig. 15). In contrast to detail D (result of the presented strategy), details A-C show significant undercuts and pronounced tool marks: signs of an unstable machining process. Recorded with Keyence VHX-7000 equipped with Dual Objective Zoom Lens VH-ZST. The vertical strips are a result of lightning and do not correspond to any surface features In conclusion, the experimental results show that the proposed machining strategy which was designed to find stable machining conditions led to a clean surface.

Discussion
Milling of Blisk blades is a challenging task due to a significant risk of chatter caused by the high dynamical compliance of the workpiece. The situation is further complicated by the influence of gradual material removal on the dynamical properties of the workpiece, necessitating modulation of machining parameters such as tilt and lead angles, stock material and spindle speeds throughout the machining process. In the specific case of Blisk blades, the geometric constraints significantly limit the possibility of avoiding chatter through the setting of tilt and lead angles. Thus, the main parameters that can be used to control the stability of the machining process in this particular setting are the stock allowance thickness and the spindle speed.
In the past, the methods of chatter avoidance in milling of thin-walled workpieces were studied extensively with the aim of designing the shape and thickness of stock allowance based on predictions of stability of the machining process to identify combinations of spindle speed, tool axis orientation and stock allowance which would guarantee a stable, chatter-free machining process. The previously published methods focused either on finding a single constant stability gap that would exist throughout the machining process, or else they would divide the machining process into multiple separate parts, each with different stock material distribution and different (constant) spindle speed. None of the previous studies, however, considered the possibility of a continuous modulation of stable spindle speed throughout the material removal process.
In this study, the final shape of the stock allowance is designed iteratively with the explicit purpose of finding a stability gap that changes continuously during the whole process of material removal. This, in turn, leads to a design of a continuous function of stable spindle speeds through which a single continuous chatter-free machining process can be achieved. As a result, the machining process which typically consists of multiple separate steps is significantly simplified.
The method presented in this study is subject to some limitations. In general, it might not always be possible to find a single stability gap that continuously changes throughout the whole machining process. In such cases, some discontinuities in spindle speed are unavoidable. But even in those cases, the presented method can be used to minimize the number of these discontinuities to simplify the machining process. It should also be noted that milling machine tools do not, so far, provide the option to continuously change spindle speed without stopping the motion axes. This functionality can, however, be realized in practice through a modification of the programming of the CNC controller unit as was the case in the experimental verification presented in this study.
Several possible extensions of the presented method could be considered in future research. In the presented method, the lead and tilt angle values are considered constant throughout the machining process. These could be added as optimization variables (under the corresponding constraints imposed by the geometry of the workpiece). Furthermore, the algorithm presented could be modified to maximize the average spindle speed, thus further increasing productivity. This could be achieved by prioritizing stability gaps at higher speeds, as well as by maximizing speeds at interpolation points in the process of spindle speed calculation. Finally, the method could be applied to other types of thin-walled workpieces, possibly leading to different designs of stock allowance distribution, such as considering the allowance thickness to be piecewise constant on a tile grid instead of horizontal strips. Any such stock allowance distribution design would need to take into account both the machining strategy and the geometry of the particular workpiece. In future research, the authors would like to focus on formulation of these alternative stock allowance distribution designs for other thin-walled workpieces (such as the blade mentioned in [34]).

Conclusions
This paper presents an innovative method of integrated optimization of stock allowance distribution and spindle speed variation in machining of thin-walled blades. The method was designed to make full use of the possibility of continuous spindle speed change throughout the finishing process.
The main results of this study can be summarized as follows: • The shape and distribution of stock allowance of a thinwalled Blisk blade can be designed in such a way that the need for discontinuous changes of stable spindle speeds is minimized throughout the finishing process. This simplifies the process and increases overall productivity. • This simplification of the machining process can be achieved through integration of stock allowance design with methods of material removal simulation and machining stability analysis in order to find a stability gap that continuously varies throughout the process of material removal. • Once the continuously varying stability gap is identified, a continuous function of stable spindle speeds can be designed which prescribes spindle speed at each point of the toolpath. This eliminates the need for multiple stages of material removal and multiple tool passes, greatly simplifying the finishing process.

Appendix A. Optimization results
This appendix contains additional information about the optimization results. Namely, the final stock allowance distribution (Fig. 19) and the stable spindle speed intervals associated with each strip (Table 5).

Fig. 19
End product of the semi-finish phase. Note the stock allowance distribution

Appendix B. Comparison of machining strategies
This appendix presents details of blade surfaces machined with the machining strategies A, B, C and D (the optimization strategy presented in this paper), Figs. 20 and 21.

Appendix C. Material removal simulation software MillVis
The material simulation software used, MillVis, uses a distance field representation of material surface. In this representation, the material surface is implicitly defined by a distance function evaluated in specific points of a uniform 3D grid consisting of cubic blocks, which are further subdivided into 3 × 3 × 3 cubic cells. Each block is assigned a flag that specifies whether it lies completely inside or outside of the workpiece, or whether it intersects the surface. The distance function values at the vertices of cells are computed only for those cells that intersect the workpiece surface. The distance to the surface at points inside such cells is then approximated using trilinear interpolation. A schematic representation of the distance field structure and the implicit surface representation is shown in Fig. 22.
During the material removal simulation process, the surface and inside blocks that are close to the tool are identified. Then, for each of their cells, it is checked whether the cell is intersected by the tool, and if so, the distance function values at the corresponding vertices are updated. In practice, the representation of the workpiece geometry by a distance field structure has been shown to be sufficiently accurate while maintaining reasonable memory requirements, as lower accuracy is usually required for larger workpieces and vice versa.
The workpiece represented by a voxel distance field can be effectively visualized using the ray tracing method, which provides high fidelity detail of the machined surface.
The definition of the tool envelope and the formulation of the cutting force is based on the framework by Engin and Altintas [58]. The ability of MillVis software to accurately predict cutting forces has been repeatedly verified by comparing its results with dynamometer data. As an example, a comparison of measured and simulated cutting forces during pocket machining is shown in Fig. 23. The corresponding experimental setup and MillVis simulation are shown in Fig. 24.

Availability of data and material Not applicable
Code availability Not applicable

Declarations
Ethics approval Not applicable

Conflict of interest
The authors declare no declare competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.