Modelling of Machining Processes

The Twin-Control product takes a holistic approach to modelling and simulation of machining processes, incorporating machine controller, machine structure, part program, part geometry, and process forces and dynamics into a single system. The focus of this chapter is on the process models used to simulate cutting forces, torques, form error, and tool dynamics throughout a part program. These predictions provide process planners with valuable information about an operation before it is executed on machine, allowing for potential issues to be identified and reduced or eliminated before the first part is produced.


Introduction
Interactions between cutting tool and workpiece are critical in any machining operation. As the tool moves through a workpiece, cutting process induces forces on both the tool and the workpiece. These forces in turn have an effect on the process and can have a detrimental effect on the machine, tool and resulting part under certain conditions. It is critical then to understand these interactions before a part is machined to avoid scrapped parts or damage to the tool or machine. This chapter first covers the development of process models used to predict cutting forces for specific part programs. The cutting force model is then used for dynamic analysis to determine the effects of tool vibration on the final part outcome.

Discrete Cutting Force Model
The 5-axis machining operations bring new challenges for predicting cutting forces, where complex tool-workpiece engagements and tool orientations make it difficult to adapt 3-axis process models for 5-axis operations. A model is developed here to predict cutting forces with arbitrary tool/workpiece engagement and tool feed direction. A discrete force model is used, in which the tool is composed of multiple cutting elements. Each element is processed to determine its effect on cutting forces, and global forces are determined by combining the effects of multiple engaged elements. The cutting force model is combined with ModuleWorks software which predicts toolworkpiece engagement regions (TWE) based on tool motion through the workpiece geometry. Cutting forces are predicted throughout complex operations by applying TWE data to the elements of the force model. The force model is validated through cutting trials in which the measured forces are compared with predicted forces during 5-axis milling.

Introduction
Cutting forces depend on the tool and workpiece material, cutting tool geometry and cutting conditions. In 5-axis milling, cutting conditions can vary considerably in process, and the varying cutting conditions can result in complex tool-workpiece engagements (TWE). Ozturk and Budak calculated TWE analytically for 5-axis ball end milling [1] and simulated cutting forces throughout a toolpath after calculating the cutting parameters at discrete intervals [2]. Although this method gives accurate results for smooth machining operations, the analytical engagement model loses accuracy when the uncut surface is more complex.
More detailed engagement calculation methods have been developed to better simulate more complex machining operations. These models operate by creating a virtual workpiece and removing any material that interferes with the geometry of a tool moved along a path. For each tool motion, the surface patches of the tool that remove material are the TWE region. In the solid model-based material removal simulation, the engagement area is derived from finding intersections between the solid models of both the tool and the workpiece [3][4][5]. Taner et al. used a planar projection strategy to determine TWE for constant feed and constant lead and tilt operations [6]. Others have represented the workpiece as a Z-map, also known as height map, a matrix/manifold of lines which are virtually cut when they interfere with the tool mesh [7]. A more advanced version of Z-map is the dexel approach [8] that can model overhangs in the geometry, thus supporting 5-axis milling. The dexel approach may be improved to so-called tri-dexel model by introducing virtual grid lines in three directions to reduce dependence on grid directionality in the geometry accuracy for any cut direction [8][9][10].
In the current model, ModuleWorks software which applies the tri-dexel model is used to determine TWE data for every cutter location (CL) point of a part program (see Fig. 4.1a). This TWE data determines which elements of a discretized tool mesh are engaged in the cut during that move. The cutting force contribution of each engaged element is then combined to determine the global cutting force for that move.

Discretized Force Model
In order to predict cutting forces for arbitrary feed direction with arbitrary TWE, a discrete cutting force model is used. The model concept is shown in Fig. 4.1b, where cutting forces on a bull nose end mill act in different directions based on the cutter position and orientation. An example of the local cutting forces is shown at one section of the cutting edge, where the local radial force F r , acts inward, normal to the cut surface, the tangent force, F t , acts in the opposite direction of the cutter motion, and the axial force, F a , acts tangent to the cut surface along to the tool profile.
The complex cut area from Fig. 4.1c is discretized into multiple elements in Fig. 4.1d. Each element has an effective cut width, b el , along the tool profile, and thickness, h el , normal to the tool profile. The global tool force is determined by combining the effects of all active cutting elements.
This section outlines the processes to determine the effects each tool element have on global cutting force.

Tool Discretization
The tool is discretized circumferentially into elements, dθ , and along the tool profile, L, into elements b el . Discretization along the tool profile and circumferentially allows the TWE of the tool to be defined by a single 2D matrix. Elements along L are created with equal length, b el , regardless of the orientation of the elements, and the element size circumferentially is dependent on the radial position of the element, r el . An example mesh for a bullnose end mill with diameter, D, corner radius, R c , and a maximum axial length of Z max , is shown in Fig. 4.2a. The mesh is created to follow the helical curve to match the shape of the cutting edge, as shown in Fig. 4.2b for zero helix angle and 30°helix tools. The tool is discretized along the tool profile, L, into N L elements, and circumferentially into N c . The mesh structure is shown in Fig. 4.2c, with L mesh indices representing concentric circles radiating from the tooltip center and extending up the side of the tool. The element cut width, b el , is the distance between two adjacent L elements. The indices indicate the circumferential position of the elements. The elements are positioned with lag angles to follow the helical curve of the cutting edge, as shown in Fig. 4.2c. By creating the mesh along the helical curve, the indices of the TWE map always correspond directly to the cutting edge, and each index corresponds to the elements of one flute at one rotational position.
Each element of the mesh has the indices, el( , L), and is defined by a set of position coordinates in Cartesian (X el , Y el , Z el ) and polar (r el , θ el , Z el ) coordinates and by an orientation angle, κ el (see Fig. 4.2a).

Tool Coordinate Systems
Three coordinate systems (CS) are used in the cutting force analysis for each element of the cutting edge; {rta}, {R T A} and {XYZ}. {XYZ} is the tool global CS, in which tool motions and cutting forces are determined. {R T A} is the tool global polar CS, describing radial, tangent, and axial directions. {rta} is the local element CS and is used to account for the orientation of the cutting elements. The rta directions correspond to the force directions associated with the cutting force coefficients (CFCs), K e,rta and K c,rta , which are fixed relative to the orientation of the cutting element, but Forces in these three coordinate systems are shown for a single element in Fig. 4.5. Ultimately, forces in the XYZ direction are required based on tool feed in the XYZ direction, but the cutting forces must first be determined in {rta} in which the CFCs are defined. Operations required to transform forces from {rta} to {XYZ}, as shown in Fig. 4.3, are discussed in this section.

Element Cutting Forces
Each tool element will contribute to the global cutting force if it is engaged in the cut. In this section, the element cutting forces are first calculated in the local {rta} CS and then transformed to the global tool {XYZ} CS.

Element Local Cutting Force
Element cutting forces are calculated in {rta} using Eq. (4.1). To determine these element forces for any feed direction, chip thickness, h el , in Eq. (4.1) is defined based on the relative feed of that element in the local rta directions, f rta (feed per tooth vector in the rta directions). By defining the element uncut chip thickness, h el , as a function of the feed vector, it is possible to calculate element forces for any feed direction from a single cutting force matrix for that element. This feature is especially convenient in 5-axis machining where the tool continually changes feed direction.
The relationship between the uncut chip thickness and the relative feed of the element, f rta , is illustrated in Fig. 4.4, where h el is equal to the distance that the element is fed in the negative r-direction, f r (the local r-direction vector always points inward, normal to the cut surface). Defining the uncut chip thickness for each element in Eq. (4.2) and combining with Eq. (4.1), the resulting expression for the cutting force for each element is given in Eq. (4.3).

Cutting Force Transformations
Cutting forces from Eq.
The second transformation accounts for the angular position of the element around the tool. The transformation matrix, T θel in Eq.  T el T θ el T κ el (4.6) Forces in {XYZ} are found by performing a coordinate transformation to the edge force vector in (4.7), and the vector field rotation on the cutting force matrix in (4.8).

Tool Cutting Forces
After transforming in Eqs.   The global cutting force on the tool for each flute position is determined using Eq. (4.11). Note if multiple flutes are engaged in the cut at once, the effects of all engaged flutes can also be combined.
Cutting forces are determined by evaluating Eq. (4.11) for each flute position. As the tool rotates through the different positions, the components of Eq. (4.11) change to reflect the engaged elements in that section.
The example in Fig. 4.7 shows the calculated cutting forces for a two-fluted, 12mm-diameter ball end mill with full radial immersion at 10,000 RPM with a cutting depth of 3 mm. For this cut, the tool is fed in the negative Y -direction at a rate of 0.1 mm per tooth, so the feed vector is constant at f XYZ {0, −0.1, 0} T . Equation (4.11) is then evaluated at each position (only positions 1 through 7 are shown), to obtain the changing cutting forces as the tool rotates through the engagement region.

Part Cutting Forces
In 5-axis machining, the orientation of the tool can change continually with respect to the workpiece. So far cutting forces have been determined in the tool CS. However, cutting forces in the part CS, XYZ P , based on tool feed in XYZ P are often more convenient as tool motions are programed in the part CS. After transforming F e,XYZ ( ) and Q XYZ ( ) in Eq. (4.12) based on the coordinate transformation from the tool CS to the part CS, T T 2P , cutting forces are obtained in the part CS using Eq. (4.13).

Cutting Trials
Two cutting tests are performed to compare simulated and measured cutting forces. The tool for both tests is a 12-mm-diameter ball end mill with two flutes and helix  . The machine used for the tests is a MAG FTV5-2500, and force measurements are collected using a Kistler 9139AA dynamometer. The CFCs used for the AL7075-T6 workpiece are found experimentally using a ball end mill mechanistic model [11], to be: K e,r 13.9, K e,t 7.1, K e,a −1.3 N/mm, and K c,r 619.9, K c,t 1014.2, and K c,a 58.2 N/mm 2 . Note that average CFC values identified experimentally are used for all elements regardless of local oblique and rake angles. A summary of CFCs for this set of trials ("M" machining tests) along with the CFCs for other trials in this chapter are provided in Table 4.1.
The first test is a stair step test where the axial depth increases from 1 to 3 mm in increments of 0.25 mm, with full radial immersion and zero lead and tilt. For this test, the spindle speed is 10,000 RPM with a feed rate of 0.1 mm per tooth in the −X P direction. The force measurements in the part CS, F XYZ,P , are shown in Fig. 4.8, along with the maximum and minimum simulated forces. A detailed comparison of the simulated and measured cutting forces is also shown in Fig. 4.8 for the minimum depth of 1 mm and maximum depth of 3 mm. It can be seen from these results that the simulated results closely match the experiment for this simple operation.
The second test is on the "M" part shown in Fig. 4.9, which is machined at 12,250 RPM using approximately 2000 CL points. During this test, the A and C axes of the machine are fixed at 21.1°and −134°, respectively, resulting in a lead and tilt   Test cut "M" character with varying cut depth of 15°when feeding in the -X P feed direction on the X p Y P plane. The rotary axes are fixed to reduce tracking errors between the simulation and the true machine motions. This part is machined over a dome-shaped base which results in varying depth over the "M" profile. Further, the fixed tool orientation results in a fixed transformation, T T 2P , but variable lead and tilt angle as the tool changes direction. Prior to the start of the simulation, the tool mesh is created, and F e,XYZ,el and Q XYZ,el are determined for each tool mesh element (these only need to be calculated once for a given tool mesh and fixed set of cutting force coefficients). Then, for each CL point, the TWE is obtained for that move using ModuleWorks Software. Figure 4. 9 shows examples of how the TWE maps obtained change throughout the operation. The TWE data for each move is then applied to obtain F e,XYZ ( ) and Q XYZ ( ), and the cutting forces are calculated using Eq. (4.11) at each angular position.
Note that the transformation between the tool CS and part CS is dependent on how the tool CS is defined in the part CS. In the step example shown in Fig. 4.8, the tool CS is the same as the part CS, so no transformation is needed (except to reverse the sign of the forces to represent forces acting on the part, where forces are calculated as acting on the tool). The transformation, T T 2P , for the 5-axis example Due to tracking errors in the tool feed velocity (despite fixing the A and C axes), the total machining time is measured as approximately 20% longer than expected. To account for this in the simulation, the feed per tooth was reduced to 0.078 mm per tooth to match the actual and simulation machine time (this effectively reduced the simulated average tool feed speed, and hence feed per tooth).
The resulting forces are plotted in Fig. 4.10a along with the maximum and minimum simulated force values at those locations in time. The results show that the simulated maximum and minimum forces closely follow the force profile throughout the operation. Further, the simulated rotation dependent forces in Fig. 4.10b also closely match the measured forces. The only places where large deviations occur are in locations where chatter occurs. While the current force model does not predict forces during chatter, this model can be used to predict whether chatter will occur. This stability prediction approach is discussed in the following section.

Cutting Force Model Summary
The cutting force model developed here was created to predict forces for complex machining operations where the tool/workpiece engagement is complex and highly variable throughout. The key feature of this model is that it treats the elements of a discretized tool as individual entities which have predetermined force characteristics (F e,XYZ,el and Q XYZ,el ) which are independent of the feed rate and feed direction of the tool. When coupled with ModuleWorks TWE software to capture effect of changing cutting conditions on the TWE, it is possible to efficiently obtain complex cutting force predictions for 5-axis milling operations.
The experiments discussed here have shown that this model is capable of predicting cutting force for a ball end mill in 3-and 5-axis operations. The stair test resulted in accurate predictions of cutting force at varying cut depths. The "M" test demonstrated that this model is able to pair with Moduleworks TWE software and accounts for a high variation of cutting conditions effectively. The predictions from this test showed close agreement between simulated and measured forces.

Stability Roadmap
Chatter is one of the major limitations to machining processes that affects part quality and productivity. It has been extensively studied [12], and different approaches have been developed to predict and mitigate chatter. Several modelling techniques such as frequency domain solution [13], semi-discretization [14] and temporal finite element [15], and time domain solution [16] have been applied to predict stability in 3-axis milling operations. These methods are used to create stability lobe diagrams (SLDs) which show the stable and unstable depths of cut at different spindle speeds.
The frequency domain solutions for flat end mills were later extended to general milling tools [17] and 5-axis ball end milling operations [18]. Additional degrees of freedoms in 5-axis milling complicate the calculation of tool and workpiece engagement which is a required input in modelling the mechanics and dynamics of the process. Ozturk and Budak [18] simulated the stability of 5-axis ball end milling operation using both single and multifrequency solution. They used an analytical cutter workpiece engagement boundary definition [1]. In some complex toolpaths in 5-axis milling, the definition of tool and process parameters such as cutting depth, step over, lead and tilt angles may not be enough to define the tool and workpiece engagements. For such cases, numerical geometric engagement calculation methods are required [19]. Ozkirimli et al. [20] calculated the SLDs by using a numerical engagement method employing the frequency domain solution.
In order to simulate stability of a complex toolpath, Budak et al. simulated stability lobe diagrams for different points along the toolpath [21] and generated 3D stability lobe diagrams. This approach is applicable for parameter selection while designing toolpath. By looking to the 3D stability lobe diagram, the most conservative cutting depth and spindle speed can be selected for the process to have a chatter-free process. On the other hand, there are practical issues with using 3D stability lobe diagrams.
In complex cutting cases, definition of cutting depth can be vague. Even if the process planner has the 3D stability lobe diagram for a given process, it will be difficult for the process planner to identify whether the process is stable at a given point. First, they need to determine the cutting depths at each point in the toolpath and spindle speed and compare them with the 3D stability lobe diagram. Moreover, once the toolpath is generated it is not possible to change the cutting depth without generating the toolpath again. For these reasons stability lobe diagrams are less practical for the visualization of the stability and changing the cutting depth without changing the toolpath. In order to visualize the stability of a given process, the stability roadmap (SRM) is proposed. It plots the stable speeds for a given part program considering the changing cutting conditions.
The purpose of the SRM is to efficiently obtain and represent stability information for complex 5-axis operations. Once a part program is created, the TWE is assumed to be fixed for each cutter location (CL) point of the program. Apart from regenerating the part program, the most effective way to reduce chatter is through spindle speed control. The SRM, thus, can be used by process planners to identify the best spindle speeds throughout the program.
The force and stability formulation required to generate SRMs are presented in this section, along with validation tests in which the model predictions are compared with measured results.

Dynamic Force Model
The cutting force matrices developed in Sect. 4.2.2 to predict cutting forces for any tool/workpiece engagement, and feed direction can also be applied to predict stability for the same operations. To simulate dynamic effects, Eq. (4.11) is modified by including dynamic effects in the tool displacement vector, XYZ . The dynamic feed vector is shown in (4.14), where f XYZ is the nominal feed per tooth vector, and XYZ and XYZ τ are tool displacements from the nominal position at the current time and at the previous tooth pass.
The cutting force matrices, Q xyz ( ), are also applied for stability prediction using the zeroth order stability model developed in [13]. For the zeroth order linear milling stability formulation, the nonlinear relation between force and feed is not included in the model. The variable component of the cutting force is expressed in the frequency domain in Eq. (4.15) [13].
The dynamic force in Eq. (4.15) has time-varying coefficients due to the tool rotation angle dependence of Q xyz ( ) (each rotation angle can be expressed as a function of time as the tool rotates). Equation (4.15) is made time-invariant by approximating the time-varying term as an average value [12]. Averaging the Q xyz ( ) over the tooth passing period, or over the range of rotation angles, , is analogous to the approximation of time-varying coefficients through the zeroth order Fourier series term in [12]. The average Q matrixes over one tool revolution, Q 0 , are calculated in Eq. (4.16), where N c is the total number of circumferential elements in the tool mesh.
After averaging to obtain Q 0 and defining XYZ (ω) as a function of the frequency response function, FRF(ω), in Eq. (4.17), the resulting time-invariant dynamic force equation is obtained as Eq. (4.18). Note here that Q xyz and Q 0 are the force matrices in the tool CS, which should correspond to the FRF measurements taken at the tooltip.
The dynamic stability is determined through the determinant of the characteristic equation in Eq. (4.19), which defines the condition at the limit of stability [12].
For calculation, Eq. (4.19) is rearranged to the form used by the eig function in MATLAB in Eq. (4.20), as done in [22].
Following the steps of reduction in [13,22], the condition of Eq.
Equation (4.21) defines the conditions when the system is at the limit of stability. Because the current system is based on a discrete model, there is no additional variable which can be solved to satisfy Eqs. (4.20) and (4.21) (such as depth of cut in [13]). As a result, stability is predicted for each new Q 0 matrix, where Q 0 changes throughout a program based on changing TWE data. For each new Q 0 matrix, the eigenvalues are found over a range of frequencies using Eq. (4.20), and the conditions of Eq. (4.22) are used to determine at which frequencies the system is stable or unstable.  Equation (4.23) is then used to determine spindle speeds at which chatter will occur, Ω c , based on the chatter frequencies. The chatter frequency is used to calculate unstable spindle speeds for multiple integers of the chatter frequency (similar to the multiple lobes of a stability lobe diagram) in Eq. (4.23). The following section describes this process in more detail.

Stability Roadmap Generation
The process used to produce a SRM from the system eigenvalues is shown for three example engagements in Fig. 4.11. For each engagement, the maximum real eigenvalues are identified over a range of frequencies, and values greater than 0.5 are considered chatter frequencies, ω c , as indicated for the third engagement in Fig. 4.11a. The chatter frequencies correspond to spindle speeds, c , which repeat for each jth lobe, according to (4.23). Here, start and end chatter frequencies, ω c, (1,2) , are identified at location where the eigenvalues cross 0.5, and these values are used to calculate start and end chatter speeds, c, (1,2) , for each jth lobe for each engagement, as shown in Eq. (4.24) and for the third engagement in Fig. 4.11b.

Stability Roadmap Trials
Stability roadmaps are generated for two example part programs and cutting trials are performed to compare the stability predictions with stability measurements. The first trial is a 5-axis ball end milling of a 3D "M" character in which the A and C axes are fixed, and the second is a continuously varying lead and tilt corner cut. The stability predictions from both examples are tested experimentally.

Stability Roadmap Trial 1
The "M" cutting trials used to test force predictions in Sect. 4.2.5 are also used here for stability prediction. The only new information required is the FRF data of the tool, which is plotted for the tool mounted on the FTV-5 in Fig. 4.12a. Figure 4.12a shows that this tool has significant dynamic asymmetry between X-and Y -directions, with significant magnitudes in the cross FRF measurements at the tooltip. As such, cross FRFs are considered in stability predictions. The resulting SRM for this operation, shown in Fig. 4.12b, is generated by following the steps shown in Fig. 4.11, using measured tool FRF data and the TWE data from each CL move. Cutting trials are run at 12,250 and 10,850 RPM. After each test, a spectrogram of the microphone signal is generated to identify tooth passing and chatter frequencies throughout the cut. Dominant frequencies, excluding tooth passing harmonics, are considered to be a result of chatter.
The results of both tests are shown in Fig. 4.13. The SRM color map indicates the predicted chatter frequency at that location, where the chatter frequency with the highest real eigenvalue, λ max , is represented for each spindle speed (there can be multiple chatter frequencies at the same speed). The results from the 12,250 RPM test (Fig. 4.13a) show that two frequencies are predicted to be dominant, at 1500 and 2300 Hz. The spectrogram results show high amplitudes at both frequencies when unstable during the first 850 CL points.  Both machined "M" parts are shown in Fig. 4.13 with locations of visible chatter on the part surface labelled. The CL points corresponding to the chatter locations are indicated both on the SRM and on the photos of the part for sections A1 and A2 at 12,250 RPM and B1 and B2 at 10,850 RPM. It can be seen that these chatter locations agree with the spectrogram measurements and SRM predictions for the first 850 points.
In both tests, the tool motion path begins at point P1 in Fig. 4.13 and follows the "M" path until it reaches point P2 (P2 corresponds the CL point 850). As the A and C axes positions are constant, the lead angle on the tool becomes negative from P2 to P3 and tip of the tool is engaged with the workpiece. One potential cause of chatter prediction errors starting at P2 is the indentation effects caused by tooltip engagement start at this point. This indentation effect, which results in underprediction of the stability, is not considered in this model.

Stability Roadmap Trial 2
The "M" stability trials presented some challenges for demonstrating the SRM concept, primarily due to tooltip engagement, and tracking errors in the feed speed of the machine tool which made it difficult to correlate simulated predictions to mea- Fig. 4.14 Edge machining trial setup with no tip engagement, and continuously varying lead, tilt, and engagement sured data. A second trial is designed in which a ball end mill machines the edge of the workpiece without engaging the tooltip, as shown Fig. 4.14a, b. As to the tool moves along the edge, the tool pivots about the tooltip, causing lead and tilt vary continuously (lead between 0°and 20°, tilt between −10°and 10°) as shown in Fig. 4.14c. The pivoting motion about the tooltip also creates continuous variation in the tool-workpiece engagement throughout the cut.
The second trial is conducted on a Starrag Ecospeed with data logging capability for the collection of true axis positions and spindle speeds throughout the cut. During the edge cutting trials, force measurements are collected using a Kistler 9139AA dynamometer. During the tests, the logged machine data and measured force data are synced so that force results can be tracked based on tool position. Further, the measured position data is applied directly to the stability model so that the stability predictions can be directly correlated to the measured force data. Note that the use of measured tool motion data does not change the SLM and is only used here to more easily compare simulation and measurement.
The CFCs used for the AL7075-T6 workpiece and the two-fluted 30°ball end mill with 6 mm radius (Sandvik R216.42-12030-AK22A H10F) are found experimentally to be: K e,r 9.3, K e,t 0.4, K e,a 0.6 N/mm, and K c,r 452.4, K c,t 955.9, and K c,a 235.7 N/mm 2 (also see "corner cut tests" in Table 4.1). Note that the CFC trials were repeated on the EcoSpeed with a new tool. Note that the tool used here is the same type but not the same tool used in the previous tests. New CFC tests are conducted with the current tool which better reflects the region of the tool used in the edge cutting trials, resulting in values which differ from the prior CFC tests.
Tool FRF measurements on the EcoSpeed showed symmetric dynamic behavior at the tooltip, with negligible cross FRFs at the tip compared with the same tool mounted to the FTV5 spindle. As such, only direct FRF measurements are used for stability prediction. However, some variation is observed in the FRF measurements before machining (but after spindle warm-up cycle) and after machining, as shown  Fig. 4.15. These measurements were taken on the same day, with the same equipment and setup, but several hours apart. The primary difference is the reduction in amplitude of the 2700 Hz mode by approximately 20% in the post machining measurement. A result of this reduction is that the SRM predictions will change depending on which FRF is used. To account for this variability, SRMs are created using both FRFs and compared with the stability trial results.
During the trials, the edge machining operation is repeated 26 times at different spindle speeds. For each test, machine axis positions, spindle speeds, and measured force data are collected throughout the operation. The measured machine data is then applied directly to the simulation model to simulate the average forces and stability at each measurement step. The measured force data is used to verify the simulated average force data, as seen in Fig. 4.16a, and to determine the regions of chatter. A spectrogram of the Z-force data is used to determine stability, where peak amplitudes at the chatter frequency are compared with peak amplitudes of the tooth passing frequency. As shown in Fig. 4.16b, the system is considered unstable whenever the chatter frequency amplitude is greater than the tooth passing frequency amplitude, and stable otherwise. Note that Z-force data is used here because it has the smallest tooth passing frequency amplitude throughout the operation, making the distinction between chatter and stability more easily recognized using the current approach.
The process for determining stability in Fig. 4.16 is repeated for all machining trials, and the results are shown in Fig. 4.17 along with the predicted SRM using both sets of FRF measurements. The predicted and measured stability results are plotted based on the position of the tool. Tool position is used here because the cycle time of each test varies due to variations in feed rate to maintain a constant feed per tooth.
The results from Fig. 4.17 show that the stability regions do not line up exactly, as there is a shift in RPM between the predicted and measured stability regions. The cause of this shift is not known for certain; however, it is possible that this shift is a result of changes in the spindle dynamic parameters as a result of spindle speed. Despite these differences and the highly transient nature of this operation, the measured regions of stability of the corner machining example closely follow  Whereas the SRM was colored to represent predicted chatter frequency in the previous example, here, the SRM colors show the predicted maximum real eigenvalues of the system, λ max (ω). The use of λ max (ω) is useful for showing the predicted severity of the chatter. This is seen in Fig. 4.18 where the experimental chatter data is plotted with a third dimension, showing the ratio of the measured maximum chatter frequency amplitude, max(Amp(ω c )), to the tooth pass frequency amplitude, max(Amp(ω t )). It can be seen that the amplitudes of the chatter vibrations qualitatively follow the form of the eigenvalues, λ max (ω).

Stability Roadmap Summary
The SRM provides an effective means of representing stability information for complex machining operations. Application of the zero-order approximation method to a discretized cutting force model allows for efficient stability prediction regardless of the engagement or the tool feed direction. When coupled with TWE simulation software, this approach can be used to represent the entire process virtually for a specific part program with a specific tool.
The experiments from this paper show that the SRM can accurately predict chatter locations, even with a tool with nonsymmetric direct FRFs and significant cross FRFs. However, the experimental results have shown that the current model is only effective when the tooltip center is not engaged in the cut.
Moving forward, the current SRM model can be expanded to include additional machining considerations, such as surface finish, surface location error, and workpiece dynamics.

Surface Location Error Model
During machining processes, cutting forces cause the tool to displace relative to the workpiece due to flexibility in the tool and/or workpiece. These relative displacements result in form errors in the part geometry, known as surface location errors (SLE). SLEs differ from chatter in that they are a result of periodic cutting forces (the forces calculated in the static cutting force module) and not regenerative effects. As such, SLEs are prevalent for both stable and unstable cutting conditions. The prediction of surface location error has been made using various techniques in prior research. In these works, the motion of the tool in response to the cutting force is predicted, and this motion is imposed onto the rotating cutting edge to predict the true surface left behind. A closed form solution to predict SLE as a function of the tool FRF and forcing function is developed in [23]. Time domain simulations have also been used to predict tool motions in more complex cases, such as run-out [24], or 2-DOF milling dynamics [25]. Others used a truncated Fourier solution to determine tool motions based on modal parameters [26]. They then simulated the full 3D "morphed" cutting edge path, which then used to model the final machined surface.
Ozturk et al. [1] developed an approach for predicting form error for 5-axis ball end milling operations. Their approach predicts form error based simulation of static tool deflections, and these predictions were shown to closely agree with single point measurements in machining trials.
In this section, the strategy used to approximate SLE based on the discrete tool model is developed. This model is able to predict 5-axis SLE for complex tool geometries, tool/workpiece engagements, based on dynamic displacement and rotation of the tool. Cutting trials are conducted to compare predicted and measure SLE for a 5-axis ball end mill operation.

Form Error Prediction
Surface location error is predicted by first simulating tool displacements based on the cutting force profile and the system dynamic parameters. These displacements are combined with the rotary motion of the tool to model the true path of the cutting edge (helical or straight). The true cutting edge locations are compared with nominal (with no relative displacement) edge locations to determine the edge position error at each location of the tool.
These steps are shown in Fig. 4.19 for an example case. In Fig. 4.19a, the simulated tool displacements are shown which is based on the cutting force for the current operation and the frequency response function (FRF) of the tooltip.
Once the tool displacements are known, they are combined with the rotary motion of the cutting edge, as shown in Fig. 4.19b. With no tool displacements, the helical cutting edges nominally follow a path that forms a cylindrical shape. The addition of the tool vibrations causes the helical cutting edges to form a new shape, which represents the actual profile of the tool. By following the shape traced by the cutting edges, which is plotted with green points, the deviation of each point can be determined and the surface location error identified at each point of the tool. In the example in Fig. 4.19b right, the tool edge trace is oriented so that the tool is feeding out of the page, and the machined surface left behind is the leftmost side of the tool profile. Nominally, the tool removes material along the cylindrical shape, leaving behind the straight, blue edge. When vibrations are included, the resulting machined surface is curved, resulting in overcut along the upper surface, and undercut at the lower surface.
In the current model, dynamic displacements are derived using a frequency domain solution. Compared with time domain simulation, the frequency domain approach is intended to be more efficient and robust. Further, this approach can be applied using system FRFs directly without the need to identify modal parameters, which is a requirement for time domain models. The frequency domain solution does not account for regenerative force effects which lead to chatter, so it is run under the assumption of stability (chatter is considered using a Stability Roadmap).
The process of determining the frequency domain solution is shown in Eqs. (4.25) through (4.27). The predicted force signal (F(t) in Fig. 4.19) is first transformed to the frequency domain using the fast Fourier transform (fft in MATLAB) as shown in (4.25).
Equation (4.26) is then used to solve for position in the frequency domain using the tool dynamic data (FRF(ω) in Fig. 4.19). The 3 × 3 FRF matrix in (4.26) allows for all direct and cross FRFs to be used when determining the tool response, although, typically only the diagonal, direct FRFs are nonzero. However, cross FRFs were used in the validation tests presented here due to their significant magnitudes.
Finally, the resulting time response of the tool is obtained in (4.27) using the inverse fast Fourier transform (IFFT in MATLAB).
The resulting time response obtained using the frequency solution is compared with time domain simulation results in Fig. 4.20. It can be seen that even though regenerative effects are ignored, the results are nearly identical using either method. The advantage of the frequency domain solution is that solutions are obtained more quickly, and determination of modal parameters it is not required. Note that this process is only valid for stable cutting conditions, and the results will differ if unstable. As such, this approach should be used in combination with the stability roadmap to ensure that the process is stable.

Surface Location Error Calculation
The tool position response is simulated over multiple cycles to increase the likelihood that the vibrations are in steady state. Once in steady state, the motions from a single tooth passing period are needed to model the tool motion. At each step of the tooth period, the instantaneous displacement of the tool is used to determine the  Fig. 4.21, starting with a single period of the tool motion. Each point in time of the tools' motion corresponds to a tool rotation position, and this displacement is applied to the cutting edge at that angle. As SLE is measured from the normal of the part surface, the SLE of each element along the cutting edge is determined by projecting the displacement vector, u XYZ , onto an edge normal vector, V n , as shown in Eq. (4.28). The resulting SLE is then recorded for that element on the cutting edge in that rotational position. This process is then repeated for each point of the tool rotation cycle. Once completed, the SLE data is represented tool mesh element using a color scale, as shown in Fig. 4.21.
SLE el( ,L) V n,el( ,L) · u xyz ( ) (4.28) The full SLE map in Fig. 4.21 shows the SLE of each element of the tool mesh; however, we are only concerned with tool elements which are both engaged in the workpiece and are also located along the section of the tool which leaves behind a machined surface. These requirements are shown in Eq. (4.29), where engagement is determined with g( , L), and "surface forming" elements are identified by projecting the feed vector, f XYZ , onto the surface normal vector. This projection gives a value proportional to the chip thickness, and it is assumed that engaged elements which have zero chip thickness are machining at locations where a machined surface is left behind. Since the tool mesh is discrete, it is unlikely to find a calculated chip thickness of exactly zero, so a minimum value is set to identify elements close to the surface formation region. A summary of the process used to identify surface forming elements for predicting SLE of the machined surface is shown in Fig. 4

Surface Location Error Trials
The SLE model is tested using a grooving operation with a ball end mill with a constant lead and tilt. The trials are run at different spindle speeds and depths of cut, and the resulting measured groove surfaces are compared with simulated surface shapes. For each set of cutting conditions, the machined surface is simulated by predicting the path of the tool cutting edges as they rotate and vibrate simultaneously, as shown in Fig. 4.23a. Note that only elements that are both engaged in the cut, and are on the edge of the tool profile when looked at in the feed direction, are used to predict the surface left behind (see elements indicated in Fig. 4.23a). Since we are comparing error results at many points along the simulated part surface, a best fit error approach is used to characterize surface errors. As the nominal shape of the ball end mill machined grooves is circular, circular fit errors are used.
Once a simulated surface is obtained, a circle is fitted to the simulated surface points (see best fit circle in Fig. 4.23b), and best fit errors are determined for each point (see in Fig. 4.23c). Even though SLE errors from the nominal surface are simulated directly, the best fit error approach allows for the simulated surface shape to be evaluated independently of a known reference (i.e., the nominal surface). This is useful for comparing simulated surface errors to the measured surface errors, which do not have an easily obtainable absolute reference. Both the simulated and experimental surfaces are evaluated using the same process based on a best fit circle.
A 12-mm-diameter ball end mill with two flutes and helix angle of 30°(Sandvik R216.42-12030-AK22A H10F) is used for this trial. The tool is mounted using a Bilz ThermoGrip T1200/HSKA63 tool holder with a tool overhang of 67.2 mm. The machine used for the tests is a MAG FTV5-2500, and force measurements are collected using a Kistler 9139AA dynamometer. The workpiece material is AL 7075 T6 and PTFE nylon. The CFCs used for the AL7075-T6 workpiece are found experimentally using a ball end mill mechanistic model [11], to be: K e,r 7.43, K e,t −2.98, K e,a −2.7 N/mm, and K c,r 128.35, K c,t 965.49, and K c,a 85.34 N/mm 2 (also see "corner cut tests" in Table 4.1). Note that average CFC values identified experimentally are used for all elements regardless of local oblique and rake angles.
During the trials, the ball end mill machines straight grooves on an aluminum and nylon workpiece. Two materials are used which have very different cutting force coefficient (CFC) values so that the resulting surface can be compared with equal cutting conditions, but different forces. Each test cut is run with a fixed feed, spindle speed, lead, tilt, and cut depth. A summary of the test conditions is shown in Table 4.2, and an image of the test block is shown in Fig. 4.24a.
After the tests, an Alicona Infinite Focus G5 is used to measure the shape of the machined grooves. This machine and an example of a surface produced during the  Fig. 4.24b, c. Surface data generated in these measurements are used to characterize the true form of the grooves. As they are all machined with a 6 mm radius ball end mill, the nominal surface generated should also have a circular form with a radius of 6 mm.
To measure the true form, virtual traces are taken on the surface data measured on the Alicona. This process is shown in Fig. 4.25, where two traces are collected, providing the true form of the groove at two locations. The trace data is then compared with a simulated surface produced by simulating the tool vibrations as it cuts.
The measured and simulated surfaces are analyzed using the same process to determine the SLE in each. First, a circle is fitted to the measured and simulated surfaces, each circle having a best fit center point and radius. In Fig. 4.26, the resulting best fit radii are shown for the simulated and measured surfaces for all 14 tests. Note  that the fitting routine only considers points along machine groove, which represent only a small portion of the full 12 mm diameter circle, causing the best fit radii to be sensitive. Despite this, the resulting best fit radii closely follow the same trend for both the simulated and measured surfaces, although the measured surfaces appear to consistently have a larger radius, approximately 50 µm greater than the simulated surface. This difference may be due to errors in the actual tool geometry, such as run-out, errors in the measurement, part vibrations, or other machine errors which are out of scope of the current simulation model. One challenge for SLE experiments is determining an absolute reference from which to measure the errors in the true surface. For example, it is difficult to define the nominal center of the tool from which to compare the measured surfaces. For this reason, the best fit circles are used as a reference to measure errors for both simulated and measured surface data.
Once best fit circles are found for the simulated and measured surfaces, the deviation of each point along the surface from its respective best fit circle is set as the best fit error for that point. In Fig. 4.27a, the points of both surfaces are plotted along with the best fit circles. The errors are recorded along the trace direction, resulting in the best fit error plots shown in Fig. 4.27a. The best fit errors can then be exaggerated to represent the shape of the machined surfaces, as shown in Fig. 4.27c. It can be seen in this exaggerated view that both the simulated and measured surfaces deviate from the best fit circle in roughly the same general form for this current example.
This process is followed for each test, and the resulting errors are shown in Fig. 4.28. It can be seen that errors in the simulated surface closely trend with the errors on the measured surface, especially in the 13,050 RPM, 0.5 mm depth tests (Tests 1-5). Chatter also had an effect on the results (chatter is not considered in the simulated surfaces). When chatter is severe (Tests 9 and 10) the surface errors deviate greatly, however, for light chatter (Tests 11 through 14), the errors still trend closely, although there is additional waviness in the measured surface.
Both aluminum and nylon are machined during the tests. The purpose of this is to compare the surfaces when there is a significant difference in the cutting forces acting on the tool. For these tests, the mean forces acting on the tool are approximately 10 times greater during the aluminum sections. As a result, the steady-state vibration amplitudes which lead to SLE are expected to reduce by a factor of 10 while cutting nylon. Due to the inability to generate full measurements for the nylon surfaces (the Alicona could not measure the full sample due to surface lighting issues), we do not  have analysis on all nylon surfaces yet. In Fig. 4.29, the partial data from the nylon samples is analyzed, and the errors are compared with the aluminum surface errors. These results indicate a significant decrease in error in the nylon surface, as would be expected with reduced cutting force. However, this result is only an indication, and data from more complete nylon measurements should be analyzed to fully validate this result.
The results from these tests show that the SLE model can provide a good indication of surface error for simple 5-axis operations using a ball end mill.

Process Model Simulation Interface
The process models discussed in this chapter have been developed in MATLAB. In order to make use of these models accessible to Twin-Control partners, a MATLAB graphical user interface (GUI) has been developed. The GUI contains all of the process models, and the interface allows users to execute any or all of the models from a single setup. This section describes the basic features of the Twin-Control process model GUI.

Process Model GUI Layout, Inputs and Options
An image of the process model GUI is shown in Fig. 4.30 with all of the key sections labelled. Upon launch, the user is first prompted to select a project folder which contains all of the process data for the current application. This project folder should contain the following items: • A database of tool geometries used (*.csv file), • A database of cutting force coefficients used (*.csv file), • An *.stl file of the stock part geometry (must be geometry before machining operation), • A part program file containing tool motion and tool number data.
Once a project folder is selected, the data is processed and the tool geometry and programmed toolpath are shown in the GUI environment. At this point, the user has several options to control the simulation. The following lists show the function of each section of the GUI as labelled in A. Users can select FRF files for up to ten tools used in the part program. The FRF data for all tools is plotted when the user selects "Show FRF Data". B. The user selects which models to run during the simulation. C. The user can control the resolution of the tool mesh and the tool path resolution at which all analysis is performed (analysis is performed at fixed distance points along the toolpath, and not an every CL point). The user can also select which section of the part program to analyze based on CL move number using Max/Min Move. D. Once options are set, the user selects "Start Analysis" to start the simulation. During simulation, the status of the simulation is updated in the "Status" box. Once the simulation is complete, the user has the option to view all tool/workpiece engagements along the toolpath by selecting "Check Engagements". E. After simulation, the results can be plotted based on the selections in the box.
Forces and torques can be plotted as tool rotation angle dependent, or as average values over one revolution. They can be plotted against simulated time or against CL move number. Finally, the type of analysis results can be selected.

Process Model GUI Outputs
When the simulation analysis is completed, the results are stored in an output file for analysis. An example of the data stored and replotted for one simulation is shown in Fig. 4.31. The average forces and torques are first shown against the simulated time. The remaining results are plotted along the toolpath so that results correspond to specific locations on the part geometry. The average force magnitudes and torques are plotted first, where the color scales correspond to the simulated values. From this example, peak force and torque loads occur at the outer corners of the "M" structure. The results for chatter and SLE are also plotted along the toolpath. These Outputs of AMRC process model GUI, showing simulated force, torque, chatter, and surface location errors graphically on the part geometry results are represented by both a color scale and dot size. For chatter results, the color indicates the predicted chatter frequency along the toolpath, and the dot size indicates the system eigenvalue along the toolpath, which correlates to the severity of chatter predicted (no dot indicates no predicted chatter). The color scale of the SLE data represents the maximum SLE value calculated within the surface generation section of the tool, and the dot size represents the range of SLE values over the surface generation section.
The method of display for the results in Fig. 4.31 is intended to allow process planners to quickly and easily interpret simulation data for a specific process. These results concisely show what issues may appear during the part program, and where they will occur on the part geometry.

Conclusions
The process models developed to model interactions between the tool, and the workpiece provide important predictions about the final outcome of an operation. The models developed here are all based on a discrete tool cutting force model which provides flexible and efficient solutions for complex machining operations. It has been shown through validation testing that these models are capable of predicting force and torque, as well as chatter and surface location error. By combining all of these separate model components into a single environment, we can efficiently simulate and view results for complex operations in a concise format. We will see in the following chapters how this system has been used to improve example machining operations from both the automotive and aerospace industries.
Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), 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 license and indicate if changes were made.
The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license 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.