<%BANNER%>

Applying Decision Analysis to Milling with System Dynamics Constraints

Permanent Link: http://ufdc.ufl.edu/UFE0041111/00001

Material Information

Title: Applying Decision Analysis to Milling with System Dynamics Constraints a New Frontier in Machining Science
Physical Description: 1 online resource (150 p.)
Language: english
Creator: Zapata-Ramos, Raul
Publisher: University of Florida
Place of Publication: Gainesville, Fla.
Publication Date: 2009

Subjects

Subjects / Keywords: bayesian, decision, dynamics, milling, optimization, stability, uncertainty, updating
Mechanical and Aerospace Engineering -- Dissertations, Academic -- UF
Genre: Mechanical Engineering thesis, Ph.D.
bibliography   ( marcgt )
theses   ( marcgt )
government publication (state, provincial, terriorial, dependent)   ( marcgt )
born-digital   ( sobekcm )
Electronic Thesis or Dissertation

Notes

Abstract: APPLYING DECISION ANALYSIS TO MILLING WITH SYSTEM DYNAMICS CONSTRAINTS: A NEW FRONTIER IN MACHINING SCIENCE For many products, milling comprises a portion of the manufacturing process, whether on the product itself, as in many aerospace structures, or one of the components used in its fabrication, such as the mold in injection molded parts. In either case, the milling process must be cost effective, which motivates process optimization. In a simplistic view of high-speed machining, a user might consider it sufficient to minimize the computer numerically controlled (CNC) tool path time. However, minimizing tool path time requires other factors that limit the material removal rate to be considered. These limiting factors include: chatter (unstable cutting), dimensional tolerances, surface finish, and tool wear or breakage. Each factor has an associated cost. If these constraints are violated, the manufacturer must either correct the mistake or scrap the part, both of which carry additional costs. Also, in some instances unwanted process behavior can cause damage to the tool, which leads to necessary tool replacement. In the end, the considerations must be combined to realize maximized profit. In this research a framework is developed to combine the limiting factors listed previously, as well as the uncertainties associated with each, into a profit optimization scheme so that informed decisions about parameter selections in the milling process can be made. This framework is based on decision of analysis, a combination of decision theory implementing Bayesian statistics and experimental analysis. Decision analysis provides methods of analyzing the milling process such that the decision maker can identify deterministic values, uncertain quantities or effects, and the aspects of the process that can be controlled such that the information can be combined appropriately. Once these considerations have been combined successfully, the effects of the user controlled quantities on the overall desirability of the process, measured in this research by profit, can be determined. In addition, the user can improve the results by performing experiments to add new knowledge of the system and diminish the effects of uncertainties. This research represents the first steps toward making this framework a reality. Initially, the milling system is characterized and organized using decision analysis and its visualization tool, the decision diagram. The effectiveness of this organization is tested using a discrete optimization on a group of test parameters. Then, treatments are applied to one process limiter, stability, in order to develop a continuous optimization and enable calculations of the value of information (maximum value a user would pay for information gain) and value of experimentation (value a user places on the information obtained from a particular experiment). Finally, a method to update information or beliefs about the system s stability condition using a Bayesian approach is detailed and tested using a numerical example. A second aspect of this research is the development of a new milling super diagram . This diagram provides a simple way to display information about many process limiters in the same parameter space. The initial diagram presents the combined information from both surface location error and stability within a user selected range of spindle speed and axial depth of cut values. This diagram displays milling process information in a format relevant to the user and tailored to the user s tolerance requirements.
General Note: In the series University of Florida Digital Collections.
General Note: Includes vita.
Bibliography: Includes bibliographical references.
Source of Description: Description based on online resource; title from PDF title page.
Source of Description: This bibliographic record is available under the Creative Commons CC0 public domain dedication. The University of Florida Libraries, as creator of this bibliographic record, has waived all rights to it worldwide under copyright law, including all related and neighboring rights, to the extent allowed by law.
Statement of Responsibility: by Raul Zapata-Ramos.
Thesis: Thesis (Ph.D.)--University of Florida, 2009.
Local: Adviser: Schmitz, Tony L.
Local: Co-adviser: Schueller, John K.

Record Information

Source Institution: UFRGP
Rights Management: Applicable rights reserved.
Classification: lcc - LD1780 2009
System ID: UFE0041111:00001

Permanent Link: http://ufdc.ufl.edu/UFE0041111/00001

Material Information

Title: Applying Decision Analysis to Milling with System Dynamics Constraints a New Frontier in Machining Science
Physical Description: 1 online resource (150 p.)
Language: english
Creator: Zapata-Ramos, Raul
Publisher: University of Florida
Place of Publication: Gainesville, Fla.
Publication Date: 2009

Subjects

Subjects / Keywords: bayesian, decision, dynamics, milling, optimization, stability, uncertainty, updating
Mechanical and Aerospace Engineering -- Dissertations, Academic -- UF
Genre: Mechanical Engineering thesis, Ph.D.
bibliography   ( marcgt )
theses   ( marcgt )
government publication (state, provincial, terriorial, dependent)   ( marcgt )
born-digital   ( sobekcm )
Electronic Thesis or Dissertation

Notes

Abstract: APPLYING DECISION ANALYSIS TO MILLING WITH SYSTEM DYNAMICS CONSTRAINTS: A NEW FRONTIER IN MACHINING SCIENCE For many products, milling comprises a portion of the manufacturing process, whether on the product itself, as in many aerospace structures, or one of the components used in its fabrication, such as the mold in injection molded parts. In either case, the milling process must be cost effective, which motivates process optimization. In a simplistic view of high-speed machining, a user might consider it sufficient to minimize the computer numerically controlled (CNC) tool path time. However, minimizing tool path time requires other factors that limit the material removal rate to be considered. These limiting factors include: chatter (unstable cutting), dimensional tolerances, surface finish, and tool wear or breakage. Each factor has an associated cost. If these constraints are violated, the manufacturer must either correct the mistake or scrap the part, both of which carry additional costs. Also, in some instances unwanted process behavior can cause damage to the tool, which leads to necessary tool replacement. In the end, the considerations must be combined to realize maximized profit. In this research a framework is developed to combine the limiting factors listed previously, as well as the uncertainties associated with each, into a profit optimization scheme so that informed decisions about parameter selections in the milling process can be made. This framework is based on decision of analysis, a combination of decision theory implementing Bayesian statistics and experimental analysis. Decision analysis provides methods of analyzing the milling process such that the decision maker can identify deterministic values, uncertain quantities or effects, and the aspects of the process that can be controlled such that the information can be combined appropriately. Once these considerations have been combined successfully, the effects of the user controlled quantities on the overall desirability of the process, measured in this research by profit, can be determined. In addition, the user can improve the results by performing experiments to add new knowledge of the system and diminish the effects of uncertainties. This research represents the first steps toward making this framework a reality. Initially, the milling system is characterized and organized using decision analysis and its visualization tool, the decision diagram. The effectiveness of this organization is tested using a discrete optimization on a group of test parameters. Then, treatments are applied to one process limiter, stability, in order to develop a continuous optimization and enable calculations of the value of information (maximum value a user would pay for information gain) and value of experimentation (value a user places on the information obtained from a particular experiment). Finally, a method to update information or beliefs about the system s stability condition using a Bayesian approach is detailed and tested using a numerical example. A second aspect of this research is the development of a new milling super diagram . This diagram provides a simple way to display information about many process limiters in the same parameter space. The initial diagram presents the combined information from both surface location error and stability within a user selected range of spindle speed and axial depth of cut values. This diagram displays milling process information in a format relevant to the user and tailored to the user s tolerance requirements.
General Note: In the series University of Florida Digital Collections.
General Note: Includes vita.
Bibliography: Includes bibliographical references.
Source of Description: Description based on online resource; title from PDF title page.
Source of Description: This bibliographic record is available under the Creative Commons CC0 public domain dedication. The University of Florida Libraries, as creator of this bibliographic record, has waived all rights to it worldwide under copyright law, including all related and neighboring rights, to the extent allowed by law.
Statement of Responsibility: by Raul Zapata-Ramos.
Thesis: Thesis (Ph.D.)--University of Florida, 2009.
Local: Adviser: Schmitz, Tony L.
Local: Co-adviser: Schueller, John K.

Record Information

Source Institution: UFRGP
Rights Management: Applicable rights reserved.
Classification: lcc - LD1780 2009
System ID: UFE0041111:00001


This item has the following downloads:


Full Text

PAGE 1

1 APPLYING DECISION ANALYSIS TO MILLING WITH SYSTEM DYNAMICS CONSTRAINTS: A NEW FRONTIER IN MACHINING SCIENCE By RA L ENRIQUE ZAPATA RAMOS A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2009

PAGE 2

2 2009 Ral Enrique Zapata Ramos

PAGE 3

3 To my family and friends, without whose love and support this work would not have been remotely possible

PAGE 4

4 ACKNOWLEDGMENTS I would like to thank all of the members of my committee for taking the time to review this dissertation and providing guidance along the way. I would like to thank my a dvisor, Dr. Schmitz, for giving me an extra push whenever I needed it and helping me when I required moments of inspiration. I would like to thank Dr. Ali Abbas from The University of Illinios at Urbana Champagne and Michael Traverso for their invaluable guidance and cooperation. I would also like to thank my lab mates at the Machine Tool Research Center for their support and en couragement at those crucial moments. Last but in no way the least, I would also like to thank my family and friends for their unwavering belief in me.

PAGE 5

5 TABLE OF CONTENTS page ACKNOWLEDGMENTS ...............................................................................................................4 LIST OF TABLES ...........................................................................................................................7 LIST OF FIGURES .........................................................................................................................8 ABSTRACT ...................................................................................................................................11 CHAPTER 1 INTRODUCTION ..................................................................................................................14 Motivation and Research Objective ........................................................................................14 Machining Science ..................................................................................................................16 Parameters and Inputs ......................................................................................................17 Measuremen ts and Mathematical Models .......................................................................18 Frequency response function ....................................................................................19 Force model ..............................................................................................................19 Stability ....................................................................................................................20 Surface location error ...............................................................................................22 Surface ro ughness ....................................................................................................23 Tool life ....................................................................................................................24 Tool path ..................................................................................................................25 Profit and cost calculations ......................................................................................25 Decision Analysis ...................................................................................................................26 Decision Tree ...................................................................................................................27 Influence or Decision Diagrams ......................................................................................28 Uncertainty Encoding ......................................................................................................30 Value of Experimentation and Value of Information ......................................................30 2 LITERATURE REVIEW .......................................................................................................32 3 APPLICATION OF DECISION ANALYSIS TO MILLING OPTIMIZATION ..................34 Decision Diagram for Milling Optimization ..........................................................................35 Initial Testing ..........................................................................................................................38 Deterministic Optimization .............................................................................................41 Uncertain Force Model Optimization ..............................................................................42 Uncertain Tool Life Optimization ...................................................................................44 Uncertain Force Model and Tool Life Optimization .......................................................46

PAGE 6

6 4 BAYESIAN UPDATING OF MILLING PROCESS LIMITERS .........................................48 Improving Knowledge about the Uncertain Milling System ..................................................48 Quantifying Stability Uncertainty ...................................................................................49 Calculating Profit at a Given Setting ...............................................................................55 Milling Optimization and Experimentation in the Presence of Stability Uncertainty ............56 Value with No Experimentation ......................................................................................56 Value with Perfect Information .......................................................................................57 Va lue of Experimentation ................................................................................................59 Bayesian Predictions of Functions ..................................................................................60 Updating Milling Stability using Brownian Distributions ..............................................63 5 SUPER DIAGRAM ................................................................................................................71 Concept ...................................................................................................................................71 Numerical Example ................................................................................................................73 Experimental Verification ......................................................................................................80 SLE Verification ..............................................................................................................82 Stability Verification .......................................................................................................87 6 CONCLUSIONS ....................................................................................................................97 Milling Optimization using Decision Analysis ......................................................................97 Super Diagram ........................................................................................................................99 APPENDIX A STABILI TY AND SLE MODEL DESCRIPTIONS AND CODE. .....................................102 Stability .................................................................................................................................102 Surface Location Error .........................................................................................................108 B DETERM INISTIC EXAMPLE CODES FOR PROFIT MAXIMIZATION .......................112 Main Code ............................................................................................................................112 Tool Life ...............................................................................................................................120 Fast Fourier Transform Code Used In SLE Codes ...............................................................122 C SUPER DIAGRAM CALCULATION CODES. .................................................................123 Data Gathering ......................................................................................................................124 Data Organization And Plotting ...........................................................................................130 D STABILITY UPDATING ....................................................................................................133 LIST OF REFERENCES .............................................................................................................144 BIOGRAPHICAL SKETCH .......................................................................................................150

PAGE 7

7 LIST OF TABLES Table page 31 Discretized machining parameters values used in the initial numerical study. .................40 32 Summary of optimization results for comparison. .............................................................46 41 Gaussian distribution parameters used in the uncertainty comparison between force model coefficients and FRF calibration coefficients. ........................................................50 42 Inputs for the Monte Carlo simulation. ..............................................................................53 43 Distribution parameters. .....................................................................................................63 44 Milling system parameters. ................................................................................................64 45 Profit calculation parameters. ............................................................................................65 51 System parameters used in the super diagram numerical examples. .............................74 52 System and operating parameters for SLE tests. ...............................................................82 53 System and operating parameters for SLE tests. ...............................................................87

PAGE 8

8 LIST OF FIGURES Figure page 11 Axial depth is labeled as b along the tool rotational axis. .................................................17 12 Spindle speed ( a). ...................................................................18 13 Feed rate ( f ) and feed per tooth ( ft). ...................................................................................18 14 Examples of end milling forces in the different coordinate frames. ..................................20 15 Sample stability lobe diagram. Stable and unstable zones are identified. .........................21 16 Under cutting SLE example in a down milling cut. ............................................................23 17 Demonstration of how surface roughness is generated in milling. ....................................24 18 Sample decision tree. .........................................................................................................27 19 Examples of the four different nodes in a decision diagram. ............................................29 110 Description of the different types of relationships that arrows express. ............................29 31 Decision diagram for the milling process. .........................................................................37 32 Part path used in the machining of the W = 100 mm cube ...............................................40 33 Stability diagram for the deterministic example with the lines of constant profit superimposed ( a = 4.5 mm and ft = 0.15 mm/tooth). .........................................................41 34 SLE contours superimposed on the stability diagram for the same case as Figure 3 3. ....42 35 Stability lobe diagram for the uncertain force model optimization at a radial depth of 3 mm and a feed per tooth of 0.15 mm/tooth. ....................................................................44 36 Profit to spindle speed relationship for a = 4.5 mm, ft = 0.15 mm/tooth and b = 2 mm ....45 37 Sensitivity to the different parameter choices with respect to spindle speed in the vicinity of the optimum for the uncertain tool life and force model exa mple ..................47 41 Comparison between (A) FRF and (B) force model coefficient uncertainties. .................51 42 A) Continuous (0, 1, 3, 3) beta distribution and B) a histogram of a randomly sampled (2000 N/mm2, 4000N/mm2, 3, 3) beta distribution for the tangential di rection force model coefficient. ......................................................................................52 43 Full set of 1000 stability lobes obtained from the random sampling of the force model coeff icients propagated through the Monte Carlo simulation. ................................54

PAGE 9

9 44 A) Histogram of the distribution of the stability limits at 18000 rpm. B) Normalized beta distribution (0, 1, 1.76, 2.56) for comparison. ...........................................................54 45 Probability of stability results using the cumulati ve distributions functions at each spindle speed. .....................................................................................................................55 46 Decision diagram and generic decision tree for milling without additional information from experimental data. .................................................................................57 47 Decision diagram and tree for optimization with perfect information. .............................57 48 Summarized decision diagram and decision tree for a single experimental updating procedure. ...........................................................................................................................59 49 Probability density functions (each spindle speed has a different PDF) for the numerical solution of the stability updating problem. .......................................................67 410 Probability tree for calculating the expected profit of a production run Pexp op, bop ; T ...............................................................................................................68 411 Probability tree for calculating the value of performing an experiment, ..........69 412 Representation of the updating of the probability of stability in a series of 12 tests after A) zero tests, B) 4 tests, C) 8 tests, and D) 12 tests. ..................................................70 51 A) Stability and B) SLE diagrams (shown in m in the contours) for a sample system. ...............................................................................................................................72 52 Sample super diagram created using the data in Figure 5 1 with an SLE user limit of 70 m. ..........................................................................................................................73 53 Effects of the users SLE limits on the feasible machining domain for errors A) greater than 30 m and B) greater than 70 m at a 0.15 mm/tooth feed rate. ...................75 54 Effects of the users SLE limits on the feasible machining domain for errors A) greater than 30 m and B) greater than 70 m at a 0.075 mm/tooth feed rate. .................76 55 SLE variation with respect to spindle speed at an axial depth of 4 mm and a feed rate of 0.15 mm/tooth. ...............................................................................................................77 56 Demonstration of the constant force phenomenon. ...........................................................79 57 SLE variation with respect to axial depth at 36000 rpm and a feed per tooth of 0.15 mm/tooth. ...........................................................................................................................79 58 Variation of SLE with respect to spindle speed at an axial depth of 4.98 micrometers and a feed per tooth of 0.15 mm/tooth. ..............................................................................80 59 Tool and holder combination used in the stability and SLE experiments. ........................81

PAGE 10

10 510 Example of changes in the spindle dynamics for the Mikron UCP Vario 600 used in this study as spindle speed increases. ................................................................................81 511 A) SLE experimental part as it would appear after test cuts, B) profile of the test squares (not to scale), and C) a picture of the finished part on the machine. ....................84 52 FRF measurements results (zero spindle speed) used as input to the surface location error model. ........................................................................................................................85 513 Measured and predicted SLE in the X direction (width) of the part as a function of spindle speed. .....................................................................................................................86 514 Measured and predicted SLE in the Y direction (width) of the part as a function of spindle speed. .....................................................................................................................86 515 Mean tool point FRF used in the stability calculations. .....................................................88 516 A) A schematic of the test piece with the wall of the test cut surface highlighted and B) a test part with several cut tests already performed. .....................................................89 517 Sample data for a stable cut at 4500 rpm and 4.5 mm axial depth. ...................................90 518 Sample data for an unstable cut at 4500 rpm and 9.5 mm axial depth of cut where the chatter frequency matches some of the content from the stable cut. .................................91 519 Sample data for an unstable cut at 4400 rpm and 9.5 mm axial depth of cut where the chatter frequency is slightly off from the stable content. ...................................................92 520 Sample data for an unstable cut at 6250 rpm and 9.5 mm axial depth of cut. ...................92 521 Sample data for an unstable cut at 7000 rpm and 9.5 mm axial depth of cut. ...................93 522 Comparison of the results of the stability exper iments with the theoretical stability lobes. ..................................................................................................................................93 523 Scanning white light interferometer image of the surface left during stable cutting at 4400 rpm and 4.5 mm axial depth. ....................................................................................95 524 Scanning white light interferometer image of the surface left during unstable cutting at 4400 rpm and 9.5 mm axial depth. .................................................................................96

PAGE 11

11 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Part ial Fulfillment of the Requirements for the Degree of Doctor of Philosophy APPLYING DECISION ANALYSIS TO MILLING WITH SYSTEM DYNAMICS CONSTRAINTS: A NEW FRONTIER IN MACHINING SCIENCE By Ral Enrique Zapata Ramos December 2009 Chair: Tony L. Schmitz Major: Mechanical Engineering For many products, milling comprises a portion of the manufacturing process, whether on the product itself, as in many aerospace structures, or one of the components used in its fabrication, such as the mold in injection molded parts. In either case, the milling process must be cost effective, which motivates process optimization. In a simplistic view of highspeed machining, a user might consider it sufficient to minimize the computer numerically controlled (CNC) tool path time. However, minimizing tool path time requires other factors that limit the material removal rate to be considered. These limiting facto rs include: chatter (unstable cutting), dimensional tolerances, surface finish, and tool wear or breakage. Each factor has an associated cost. If these constraints are violated, the manufacturer must either correct the mistake or scrap the part, both of which carry additional costs. Also, in some instances unwanted process behavior can cause damage to the tool, which leads to necessary tool replacement. In the end, the considerations must be combined to realize maximized profit. In this research a framework is developed to combine the limiting factors listed previously, as well as the uncertainties associated with each, into a profit optimization scheme so that informed decisions about parameter selections in the milling process can be made. This

PAGE 12

12 framework i s based on decision of analysis, a combination of decision theory implementing Bayesian statistics and experimental analysis. Decision analysis provides methods of analyzing the milling process such that the decision maker can identify deterministic values, uncertain quantities or effects, and the aspects of the process that can be controlled such that the information can be combined appropriately. Once these considerations have been combined successfully, the effects of the user controlled quantities on th e overall desirability of the process, measured in this research by profit, can be determined. In addition, the user can improve the results by performing experiments to add new knowledge of the system and diminish the effects of uncertainties. This resear ch represents the first steps toward making this framework a reality. Initially, the milling system is characterized and organized using decision analysis and its visualization tool, the decision diagram. The effectiveness of this organization is tested us ing a discrete optimization on a group of test parameters. Then, treatments are applied to one process limiter, stability, in order to develop a continuous optimization and enable calculations of the value of information (maximum value a user would pay for information gain) and value of experimentation (value a user places on the information obtained from a particular experiment). Finally, a method to update information or beliefs about the systems stability condition using a Bayesian approach is detailed and tested using a numerical example. A second aspect of this research is the development of a new milling super diagram This diagram provides a simple way to display information about many process limiters in the same parameter space. The initial diagram presents the combined information from both surface location error and stability within a user selected range of spindle speed and axial depth of cut

PAGE 13

13 values. This diagram displays milling process information in a format relevant to the user and tailored to the users tolerance requirements.

PAGE 14

14 CHAPTER 1 INTRODUCTION Motivation and Research Objective The widespread use of highspeed machining in recent decades has led to a significant body of research on issues that limit its productivity. In traditional m achining, conservative machining parameters are often applied. In this case, it is often possible to control the heat produced by the cutting process with coolant /lubrication and tool wear may be kept to a manageable level. Additionally, process stability is considered to be limited by a single width of chip and vibrations during cutting were reasonable unless chatter occurs or the selected feed rate is too high. However, with the increased capabilities of modern machining centers and tool materials/coating prior simplistic machining parameter selection strategies do not typically yield optimized productivity. With increased spindle speeds and feed rates, improved process performance (both stability and part accuracy) can be achieved by an understanding of the underlying process dynamics. Also, increased speeds lead to increased energy and heat generation which can cause rapid tool wear, effects which can be relieved with process monitoring and coolant/ lubrication. Finally, new production strategies (such as thin wall, monolithic machining to replace sheet metal build ups and finish machining of near net shape forgings) requires that the final stress state and subsurface damage of the workpiece be considered in parameter selection. These topics have received significant individual attention recent in research efforts [1 3]. However, it remains mostly the users responsibility to obtain all available information and evaluate it to plan their machining process. To aid in this difficult task, researchers have al so undertaken efforts to combine different aspects of machining and perform optimizations considering various limiting factors. In metal cutting, typical optimization objectives include

PAGE 15

15 minimizing cost or maximizing material removal rate. Typical constrain ts have included tool wear and machine or spindle power. However, these optimizations typically exclude the effects of machining dynamics and its effects on process stability and accuracy, both of which affect cost and material removal rate. Chatter (unsta ble cutting) leaves a rough undesirable surface which is typically out of tolerance ; further, the large vibration and force magnitudes can lead to tool or machine damage Even during stable cutting, the forced vibrations due to the tool periodically impact ing the workpiece can cause significant inaccuracy in the location of the machined surface. The goal of this research is to provide a framework to simultaneously consider these factors when selecting an optimized parameter set for milling using profit as the objective function. The users knowledge or state of information, and uncertainties also need to be considered during optimization. This research effort focuses on incorporating this information about the milling process into a comprehensive comparison and optimization scheme. The approach enables the user to propagate uncertainties in inputs and process limiters through to the profit comparison. Even though methodologies for propagating uncertainties from the inputs to the output (comparison quantity) are available, this propagation has not been incorporated into milling optimization. This research addresses these issues through the combination of two disciplines: machining science and decision analysis. Machining science identifies the required inputs, process dynamics and limitations, uncertainties, and mathematical relationships between important quantities and effects. Decision analysis provides the methods to relate the machining science data and quantify and propagate uncertainties through models t o the comparison quantity. In addition, decision analysis provides methods for improving the accuracy of the state of

PAGE 16

16 information (decreasing uncertainties) by including new information obtained from experiments or models and verifying the optimality of the final results. The following sections supply background information on machining science and decision analysis, including basic concepts, key variables and quantities, properties, and algorithms involved in this framework and analysis. Machining Science The milling parameters, properties and mathematical methods used in this research are described here. There are two key questions machinists/process planners must answer when evaluating and performing any milling operation: 1) Will the cut be stable? and 2) Will the resulting part dimensions be within acceptable tolerances? The stability of a cut depends on the users choice of milling conditions, the system dynamics, tool selection, and part material and geometry. Satisfying the required part toleranc es depends on both the ability to follow the commanded contour and the effects of the forced vibrations during cutting. The contour accuracy depends on part path planning, computer numerical control (CNC) capabilities, machine thermal errors, and tool sele ction and measurement. These aspects of machining accuracy have been studied extensively and, as a result, these types of errors are generally kept to acceptable limits in modern commercial machining equipment [4]. However, the forced vibrations of the tool and/or part during the machining process also affect the accuracy of the machined surface and can vary significantly from one application to another, even on a single machining center. In traditional low speed machining processes, the machining dynamics are not considered a significant issue because spindle speeds and depths of cut are conservatively selected. Therefore, the system vibrations are typically small. However, with the increased use of highspeed machining (with spindle speeds of tens to hundreds of thousands of rpm in some cases and associated allowable increased in axial depth of cut) and the increasing requirements for smaller parts with tighter tolerances, errors caused by machining dynamics can become significant. In [4]

PAGE 17

17 it was shown that error levels related to machine accuracy were one to two orders of magnitude less than the errors caused by the process dynamics for an example commercial machining center. Subsequently, this study focuses on the errors and limitations caused by system dynamics. The following subsections describe the milling process and models relevant to the new methodology, as well as all the inputs used to define the process and models. Parameters and I nputs The behavior of a dynamic system depends on the force input and its dynamic characteristics. In milling, the force depends on a combination of user defined parameters and the interaction between the tool and the workpiece. The user defined machining parameters are the axial depth of cut, the radial depth of cut, the feed rate, and the spindle speed. In addition, the user (machinist) must select the cutting tool. Tool selection affects the machining process because material composition, tool geometry, and coatings all influence the forces that are developed between the tool and workpiece. Definiti ons of these parameters are listed next. The axial depth of cut is the length of the tool parallel to the rotational axis of the spindle engaged in the cutting operation ( b in Figure 11). In order to avoid damaging the tool and system, this parameter sho uld never exceed the flute (cutting edge) length of the tool The helix angle of the tool ( ) is the angle the flutes make with the rotational axis (Figure 11). Figure 11. Axial depth is labeled as b along the tool rotational axis. The helix angle is also show in the figure. Radial depth of cut ( a in Figure 12), is the portion of the tool diameter (d) that is engaged in the cut; it cannot exceed the tool diame ter. This dimension is defined parallel to the feed direction and perpendicular to the spindle rotational axis. In addition, this parameter can be expressed as a percentage of the diameter engaged in the cut, also called the radial immersion. For example, the radial immersion is 100% in a slotting operation, where a = d. b

PAGE 18

18 Spindle speed (also shown in Figure 12 as usually expressed in revolutions per minute or rpm. Figure 12. Spindle speed ( a); the rotational axis of the tool is into the page. The feed motion of the tool can be described by two related parameters: the feed per tooth and the feed rate. These parameters are related using equation 1 1, = (1 1) where ft is the feed per tooth in units of length per tooth engageme nt, f is the feed rate in units of length per time N is the number of flutes (or teeth) on the cutter, and is the spindle speed in rpm. Feed per tooth is the distance the center of the cutter advances between tooth engagements and is not necessarily the thickness of the chip. Figure 13. Feed rate ( f ) and feed per tooth ( ft). Measurements and Mathematical Models Users need to relate the inputs and machining dynamics to the phenomena that are observed when milling. Measurements and models can satisfy this requirement. In order to better understand the milling system, information can be collected and analyzed. First, the systems dynamic response, described using the frequency response function, is measured or modeled Second, this information is related to the observed phenomena using predictive models. These ft f a

PAGE 19

19 models accept the inputs and use this information to describe the system behavior. For milling, relevant models are the force model, the stability model, and the surface location error mod el. Frequency response function In any dynamic system the time dependent motion, or vibration, is determined by the systems dynamic response and the force input and initial conditions. In milling, the system dynamics include the combined machine tool holder workpiece system. The system dynamics are generally expressed as the complex valued ratio of the displacement to the input force over a selected frequency range; this ratio is the frequency response function (FRF). The FRF can be modeled or measured ; a combination of the two is also possible [ 5]. The FRF is determined at the point of interest (the free end of the cutting tool, which usually represents the most flexible point in the system). Measurements are often obtained using an instrumented hammer to excite the structure and a vibration sensor, such as a vibrometer or accelerometer, to record the resulting motion. Force model An accurate force model is critical because the cutting force directly influences the system vibrations, which define the process stability and part errors obtained from forced vibration, or surface location errors. In this work, a force model is applied t hat uses a combination of chip area and edge effects as shown in equations 1 2 and 13, where, for each tooth in contact with the workpiece, Ft is the tangential force component, Fr is the radial force component, ktc is the tangential cutting force coeffic ient, krc is the radial cutting force coefficient, kte is the tangential edge coefficient, kre is the radial edge coefficient, and h is the instantaneous chip thickness [ 6]. See Figure 1 4. Note that, because the tangential and radial direction forces chan ge as the tool rotates, it is often convenient to project the forces into the fixed x and y directions, where x corresponds to the feed direction. The value of h is a function of cutter angle, and can be

PAGE 20

20 approximated as = sin ( ) The cutting force coefficients in milling are specific to the tool geometryworkpiece material and can also be a function of spindle speed and commanded chip thickness. F orce coefficients may be obtained from several different data sources, including transformation of cutt ing geometric and physical parameters (shear angle, friction, and forces measured during orthogonal cutting, for example), tabulated data, and mechanistic identification, where the forces are measured and the coefficients are determined using a linear regression over multiple feed per tooth values [6] For this research, the mechanistic approach is applied. = + ( 12) = + ( 13) Figure 14. Examples of end milling forces in the different coordinate frames. Up milling and down milling cases are also identified. Stability In milling, process stability depends on the phasing of vibrations on each successive flute engagement. As the cutter engages the workpiece, forces are generated which cause vibrations. These vibrations, in turn, are imprinted on the workpiece surface. If the vibrations are in phase with (i.e., match up to) the surface left by the previous tooth engagement, then the chip thickness remains rela tively the same. When the vibrations satisfy this in phase condition, the forces and vibrations in the system persist in a steady state or stable condition. On the other hand, if the vibrations from one tooth engagement to the next do not match up or are out of phase the chip thickness varies and the force amplitude varies. The force variation excites the system dynamics F t F r Down milling Up milling F x F y x y

PAGE 21

21 and causes vibrations which lead to subsequent changes in the chip thickness. This feedback mechanism can lead to unstable behavior, exhibited as self excited vibrations or chatter [ 711]. This system is described using a set of time delayed differential equations (one for each force direction in consideration) where the time between tooth engagements is the amount of the delay. T he pha se difference between engagements depends on the spindle speed and axial depth of cut acts as the feedback system gain. Figure 15. Sample stability lobe diagram. Stable and unstable zones are identified. Stability information can be presented using different input parameters (axial depth, spindle speed and radial depth). Generally, the user selects a particular radial depth at which the stability calculations are to be completed and then a stability lobe boundary is predicted in the spindle speed axial depth domain. The graphical representation of this relationship is a boundary between stable and unstable cutting zones (displayed as a function of spindle speed and axial depth) It is referred to as a stabi lity lobe diagram (Figure 1 5). Any operating point above the boundary is considered unstable, while all the points below are considered stable. Stability can be predicted in a number of ways, including temporal finite element methods [ 12], time domain si mulations, and frequency domain methods [1315]. For this research, the frequency domain approach is applied because it enable s the use of FRF measurements (or predictions) directly, without the need for modal fitting [16]. This decreases both the processi ng Spindle speed Axial depth Unstable (chatter) Stable

PAGE 22

22 time and the cost of modeling. Modal fitting usually needs some adjustment by the user in order to match measurements, which requires time, experience, and knowledge of the dynamics and mathematics involved in the modeling. In addition, the frequency domain method is analytical and, therefore less computationally expensive. In this research, the stability diagram is considered both as a deterministic discrete property and a continuous uncertain quantity varying according to available information and the user specified risk tolerance. Surface l ocation e rror Milling is an interrupted cutting process in which the teeth on the cutting tool repeatedly engage the workpiece and cut away small chips. Therefore, even under stable cutting conditions (which do not exhibit self excited vibrations) the tool experiences periodic forced vibrations. The magnitude and phase of these vibrations depend on the process parameters. These parameters include the dynamic response of the system (represented by the FRF), the exci tation frequency (which depends on the spindle speed and the number of teeth on the cutter), the cutting force model coefficients (dependent on the tool and workpiece combination), the radial and axial depths of cut, the feed per tooth, and the cutter s he lix angle. The location of the individual cutting edges as they enter (up milling) or exit (down milling) the cut as the tool vibrates determines the final location of the machined surface. The difference between the location of the machined surface and the commanded location is called surface location error (SLE). Depending on the excitation frequency and its relationship to the system natural frequencies, the tool may remove less material than commanded (undercutting) or more material than commanded (over cutting). Figure 16 shows an example of undercutting SLE in a climb (down) milling cut Like stability, SLE predictions can also be obtained from a number of approaches, including temporal finite element methods [1 7, 18], time domain simulation, analytical nonlinear

PAGE 23

23 dynamics, and frequency domain methods [19]. In order to use the same inputs as the stability calculations and decrease computational burden, the frequency domain solution presented in [19] is used throughout this research effort. SLE may be assumed deterministic, given tool selection and machining parameter selection or the effects of input uncertaint ies may also be considered. Figure 16. Undercutting SLE example in a down milling cut. Surface r oughness During machining, the teeth on the cutter enter and exit the cut while rotating about the spindle axis and being fed into the material. Each tooth on the cutter therefore follows a trochoidal path or more specifically, a prolate cycloid where the circumference of the roll ing circle is equal to the feed per revolution and the extended radius is equal to the radius of the tool The parametric set of equations used to describe this path is: = 2 sin ( 14) = 2 cos ( 15) where the radius of the rolling circle is = 2 is the angular position of the cutter tooth, and N is the number of teeth or flutes on the cutter. However, because the feed per tooth is very small relative to the tool radius, it is generally sufficient to approximate the trochoidal path using two offset circles (as approaches zero, the x and y equations approach the parametric equations for a circle) to represent the variation in chip thickness with cutter rotation. Since these circles are offset, the surface left behind is not perfectly flat; there are small peaks (cusps) left behind as shown in Figure 17. Actual Commanded x y

PAGE 24

24 Various statistics, such as the average roughness Ra, are used to describe the smoothness of the machined surface. Surface roughness may be approximated using equation 1 6 where the plus sign is used for up mill ing and the minus sign is used in down milling [ 20]. = 232 2 ( 16) Figure 17. Demonstration of how surface roughness is generated in milling. Note that the feed between tooth engagements is exaggerated in the figure to display the concept clearly. Tool l ife The milling process creates two distinct types of interactions between the tool and the workpiece: impact as the tool enters the cut (the high energy impact can catastrophically damage the tool teeth) and shearing as the chip is cut and removed from the r est of the material (the main cause of progressive tool wear over time). Repeated sliding of the chip over the rake (cutting) face of the tool at high pressures and temperatures can cause a crater to form (crater wear) and rubbing of the relief (flank) edge of the teeth with the cut material causes the tool edges to wear (flank and notch wear) [ 21]. A traditional approach to quantifying tool wear involves selecting one of the wear manifestations and applying a preselected limit (for example, a 0.5 mm flank wear width maximum). Then, nominally identical tools are used under different operating conditions until the preselected wear limit is exceeded. The time required to reach the limit is referred to as the tool life. The test data is then used in an empirical (Taylor type [ 22]) model to Axial direction Feed direction

PAGE 25

25 relate tool life to the operating conditions; see equation 1 7 as an example where T represents the tool life in minutes, v is the surface speed of the cutter teeth (often used due to the dependence of cutting temperature and wear to this quantity), and ft is the feed per tooth. The parameters used to fit the data to the model are the exponents p and q and the constant C Additional parameters can be added to the model as required by the user. = where = 60 d is the tool diameter, and ed in rpm ( 17) In this research tool life is treated as an uncertain quantity and predictions are expressed as a distribution. This is due to the fact that, even under highly controlled experimental co nditions, the measured tool life may vary for any given operating condition. Tool p ath The tool path that the machinist defines (typically using computer aided design (CAD)/ computer aided manufacturing (CAM) software) is a deterministic quantity essentia l to calculating the costs of machining. The CAD/CAM software uses a solid model of the desired part to prescribe the tool motions required to create the final product from a selected stock geometry. Once the path of the tool is defined and the operating c onditions (spindle speed, feed rate, axial depth, and radial depth) are selected, the time required to machine the part is set. The total operating time per part ( tm or machining time) is composed of both cutting and noncutting times. Only the cutting tim e ( tc) is considered for tool wear calculations. Profit and c ost c alculations The profit obtained from any given machining operation is the difference between the money required to perform the operation (cost) and the money received from the sale of the pa rt (revenue). Revenue is generally user defined. The company selling the product performs market research, determines the value of the product, and decides on a market price. The cost to produce

PAGE 26

26 the part, on the other hand, includes many factors beyond the control of the company, such as the cost of materials and equipment, arrangement and storage of materials and equipment, transportation, taxes, payment and training of employees, design costs, maintenance, and other fixed costs. Costs specific to the mach ining process include those associated with setup time, the machine operation (power, lubrication, cooling supplies, climate control, etc.), disposal of used metal and coolant, tooling, and others [23]. For this research, the cost equation described in equ ation 1 8 is applied [ 21], where rm is the cost to operate the machine per unit time, tch is the time to change a tool, and Ct is the cost of a tool; tm and tc are determined from the final part path and depend on the machining parameter settings. Note that the equation represents the cost to produce a single part. Using the same cost rate for a tool change as for machine operating time is a simplification that depends on how the value for tool changing time is selected. For example, a single machinist may be responsible for changing worn tools on several machines requiring a lower cost per part than having a machinist at each machine ( ) = + ( + ) ( 18) Decision Analysis In a deterministic model, the best result is only judged by how much improvement it provides over other options. However, no process is fully deterministic. F or example, the microstructure of two materials can vary sufficiently to give them slightly different mechanical properties. In milling, this c ould translate to different cutting coefficients This variation could make the difference between the best possible point and a point that violates constraints. These uncertainties lead to the use of statistics in modeling and prediction efforts. Here, a combination

PAGE 27

27 of decision theory (which implements Bayesian statistics) and experimental analysis is applied. This combination is known as decision analysis. Decision analysis provides methods to organize informa tion, identify which aspects of the process may affect others, visualize the entire process, separate the process into its most essential parts, enable decision making by providing results and risk factors associated with the result, and enable updates to decisions given new information. The model(s) used in the process are developed in collaboration between the user (the person with process knowledge) and the analyst (the person acquainted with statistical modeling). This can lead to simpler, more effectiv e models. Decision T ree Figure 18. Sample decision tree. The decision tree is a basic tool in decision analysis. It allows the user to identify outcomes of decisions, including the effects of uncertainties. In a decision tree, each branch is identified with either the decision that was made or the result of the uncertain process as well as the probability of that part icular event as a result of the decision or uncertainty that led to it. Decisions are typically represented as rectangular nodes, while uncertainties are identified as Result 6 Result 1 Result 4 Result 2 Result 3 Result 5 U 1 U2 D P1 P2 P3 P4 P5 P6 P7 P8

PAGE 28

28 oval nodes. An example is presented in Figure 1 8, where D denotes a decision, U1 and U2 are uncertainties, and P1 through P8 denote different probabilities of occurrence in each branch. Note that the number of outputs from decision and uncertainty nodes is set to the number of possibilities that the user decides to consider. Therefore, the t ree can be very large depending on the process. Influence or D ecision D iagrams While a decision tree presents a wealth of information, it can sometimes become too complicated. One example is a process that includes recursive statements (similar to for and while loops in computer programming). In instances such as these, a decision diagram offers an alternative tool. A decision diagram is a graph based representation of all the processes, sources of information, decisions, and calculations as a unidirectional flow of data [2 4]. The diagram is composed of four nodes types connected by arrows. The four node types are decision nodes, uncertainty nodes, value nodes and deterministic nodes (Figure 19). Decision nodes represent a set of available choices that the user makes; these choices may be in the form of discrete sets or continuous variables. They are represented as rectangles in the diagram. Uncertainty nodes, represented by ovals, summarize a discrete or continuous variable that cannot be known or controlled perfectly. Double ovals are deterministic processes or functions that present a direct result of the inputs. A value node represent the final objective of the decision being made; it is the value by which all alternatives are compared. Typically, there is a function that : 1) relates all the different processes and conditions to be satisfied in the analysis ; and 2) generates a relevant comparative quantity. A value node is generally displayed as a polygon, such as a diamond, hexagon or octagon.

PAGE 29

29 Figure 19. Examples of the four different nodes in a decision diagram. The arrows that connect these nodes have different names depending on the type of relationship they express (Figure 1 10) Relevance arrows flow fro m one uncertainty node to another and indicate that the distribution of the receiving node may be conditional to the source uncertainty node. Information arrows flow from decisions or uncertainties into a decision to provide the user with more knowledge (t o make the necessary choice). If the information comes from an uncertainty, then that uncertainty is known prior to making the decision. If the information comes from another decision, it means that the information from the previous decision should have so me weight in the choice the user is about to make. A functional arrow indicates the inputs that flow from a decision or uncertainty into a deterministic node and that the output depends directly on those inputs. Finally, the arrows flowing into the value node are referred to as value arrows. Figure 110. Description of the different types of relationships that arrows express. A B C D 2 D 1 E D D C Informational arrow Relevance arrow Functional arrow Influence arrow F Uncer tainty Decision Value Deterministic

PAGE 30

30 Uncertainty E ncoding Uncertainty encoding is used to assess a users initial belief and risk tolerance regarding a random variable for which uncertainty is known to exist, but the limits and/ or distribution of the uncertainty are not known. The data used for this evaluation can have many sources, including user knowledge / experience and previously recorded data. There may be issues such a s user bias, but these can be adjusted and updated as more data about the system is collected. For this research, classic probability encoding is used to determine the initial settings for a three degree discrete random variable during the initial modelin g efforts. In this procedure, three values of the random variable are chosen according to the current level of information and the users belief such that they represent the low, medium and high levels of the distribution [ 2530]. These are assigned at the 90%, 50% and 10% fractiles of the assumed distribution according to the cumulative probability distribution. These points are next assigned to a discrete probability mass function such that they have 25%, 50% and 25% probability of occurrence. These value s are then used in the discrete probability choices. For completely continuous models, only the initial probability distribution assessment is required and the discretization is unnecessary. Value of E xperimentation and V alue of I nformation Value of inform ation and value of experimentation are two quantities used in decision analysis to determine the desirability of the results of a set of inputs and decisions. The value of information denotes the maximum value that a user would pay to obtain perfect inform ation about a particular uncertainty. Perfect information means that there is no uncertainty in the variable; this quantity is usually likened to asking a clairvoyant (all knowing) entity for information. Reducing uncertainty to zero is impossible in pract ice so this ideal quantity is used as a reference point by which to measure the information gained through other sources, such as modeling and experimentation. Information obtained from other sources is imperfect. Therefore,

PAGE 31

31 its value (the amount the user is willing to pay to obtain that information) is always less than that of perfect information. A comparison of the two values can then be used as an optimization criterion.

PAGE 32

32 CHAPTER 2 LITERATURE REVIEW This literature review presents previous research in the area of milling optimization. The objective is to present different approaches to the problem in order to compare and contrast process and o ptimization limiters, objective functions, and optimization routines. Machining optimization has been pursued by researchers for well over 100 years [ 31]. The thrusts of the research efforts are varied, but can be loosely categori z ed by the objective func tions including costs of production [ 3243], rate of production [ 31,44 47], and profit [34,47]. Some of the earliest examples of research on machining optimization are efforts by Taylor in 1893 and 1906 [ 22,48] where the economic ramifications of tool wear and machining rate in turning we re considered. As knowledge of the machining processes increased, more detailed optimization routines emerged. In the 1950s 1970s, machining optimization research and practice focused on minimizing cost and maximizing production rates as a function of surface speed and its exponent from the tool life model [49,50] in single pass machining, and optimal tool life, number of passes, and failure costs for multi pass machining [ 51]. As this research area progressed, more aspects and variables, such as the cutting parameters [ 52] and sensitivity to these parameters [5 3], were considered in the analysis until in 1964 the use of computers to select an optimum was suggested [5 4,55]. The power of computers has increased significantly over the years, allowing for increased complexity in the optimization analyses. This has enabled statistical studies of the process to increase in number and complexity. For example, in 1966 several papers were written on machining optimization (by a selec t group of authors) which explored different facets of the problem including the linearization of the Taylor tool life equation [56], the use of profit as an optimization criteria in machining [5 7], and the use of Bayesian statistics for the machining

PAGE 33

33 opti mization and the design of experiments [58 ]. Also with the advent of computer numerically controlled (CNC), automated part path planning became a significant aspect of optimizing the machining processes [5 961]. Most of these studies concentrated mostly on turning machining operations and many more studies abound (see [ 31] and [ 62] for more details), but studies on milling optimization are more scarce, especially for highspeed milling. More recent studies focused on milling have taken advantage of increased computational power to search larger parameter spaces using global optimization techniques. Some focus on specific aspects of the milling process such as surface roughness [60], tool life/availability [ 38, 64], and part path planning [5961], was appli ed. However, researchers have begun to look at milling holistically and are considering multiple process limiters in their algorithms. Some approaches include machine limitations [4 6, 65, 66], while others also include some component limitations (surface r oughness and burrs for example) with the machine limitations [4 4,6769]. Hybrid experimental and theoretical approaches are also used such as the Taguchi method and principal component analysis combination [ 70]. However, very few researchers have considered process dynamics limitations including stability and surface location error in milling optimization schemes One example is the multiple objective optimization scheme presented in [ 71]. This research presents a new framework in which machine limitations, part constraints and the limitations presented by process dynamics can be combined to define a clearer global picture of the machining process.

PAGE 34

34 CHAPTER 3 APPLICATION OF D ECISION A NALYSIS TO M ILLING O PTIMIZATION The milling process is used in manufacturi ng facilities from the small job shop to large aerospace companies. Within the industry, there are generally two main types of users: traditional machinists and academic machinists. The traditional machinist uses conventional methods, machining parameter c ombinations, and tools identified by experience. Typically the uncertainty in obtaining acceptable machining behavior is low because this group tends to select conservative machining rates and parameters. They understand that they may not be using the mos t efficient parameters possible, but they know that their choices usually produce acceptable results and, if not, they adjust them using trial and error. Their information comes from actual machining experience. The academic machinist relies on some level of experience, but also applies models and optimizations to upgrade the performance of the milling operation. When the models, predictions and optimizations are correct, the performance can be greatly improved. However, models naturally include some degre e of uncertainty in the predictions. The academic machinist understand s this, but is willing to accept the risk and adjust the process, if necessary, to remedy any problems. The requirement for process changes is typically based on machining performance re quirements and modifications are made based on experimental results. The traditional and academic methods have both pros and cons and may even contradict one another. Therefore, some combination of the two approaches might prove useful. Decision analysis can link the two sources of information for milling optimization. The idea is to combine the predictive information obtained from the models with the knowledge and experience from the machinist and data from the process. For example, t he traditional machini st can provide the starting point and present the model with useful initial information, such as a reasonable choice

PAGE 35

35 of tools for a particular material, an initial operating point and accuracy information when machining using those parameters. The academic method would then apply that information and add to it. New information could include measurements of the system dynamics and models of the process dynamics, for example Ideally, the combined data would allow the machinists to select operating parameters that minimizes cost according the user s or customer s needs. Decision analysis organizes information and combines it considering different levels of uncertainty and risk tolerance using one value, such as material removal rate, production cost, or expected profit, for comparison of the collection of decisions and resulting values from each combination. Decision analysis enables the machinist to combine process calculations with limitations, including stability, surface roughness and surfac e location error. This capability differentiates the research described here from previous efforts (as presented in the literature review). The first step in the decision analysis process is to logically combine the multiple sources of information. As previously described the decision diagram can be used to visualize the entire process and complete this activity Decision D iagram for M illing O ptimization The decision diagram used in this description of the milling process contains ten nodes. There are two decision nodes, one for tool selection and one for the machining parameters. The tool selection node includes all aspects of the tool geometry, material composition, coating, and the frequency response function ( FRF ) of t he tool when clamped in the holder and spindle (tool holder and spindle information are also included). Previous research [ 72] shows that there is typically low uncertainty in measured FRFs. Therefore, in this research the effects of the FRF measurement un certainty are neglected However, if necessary, the FRF can be included as a separate uncertainty node. The machining parameter decisions include the radial depth of cut a, the axial depth of cut b, the spindle speed ft.

PAGE 36

36 There ar e two main sources of uncertainty considered in this model. The first is tool life. The Taylor type tool life model is based on a fit to data from multiple wear tests while varying the machining parameters. Acceptable accuracy within a large domain requires many tests. However, even for tools that are nominal ly identical, it has been shown that their rate of wear may differ due to imperfections in the tool manufacturing process or simply because the workpiece material, though assumed homogeneous, is not. The other primary uncertainty source is the force model coefficients used in the process dynamics calculations. Although there are several sources from which force model coefficients can be estimated, in this research the mechanistic approach is applied. As described in [ 6], cutting tests are used as the data source for a least squares fit (linear regression) to obtain the necessary coefficients. As with tool life, material and cutter combinations that should be identical provide slightly different results wh ich are presented as uncertainty in the coefficients. Five deterministic nodes are presented in this analysis. Two of these nodes are user defined, but not by the machinist and, therefore, they have no inputs. The first is labeled as design and contains a ll the information about the part to be machined: dimensions, tolerances, required surface finish, and material composition. The second is the revenue from the sale of the part. Sources for this information are varied (consumer and market studies for exam ple), but are not modified by the machinists choices. The last three deterministic nodes all involve calculations and measurements. The information in the tool path node can come from a calculation by the user; however, most parts are sufficiently complex so that the tool path generation and calculations are left to CAD/CAM software. No uncertainty is considered in these inputs since t he software calculates path length and cycle times. The user may extract or recalculate the se times (according to part path s and feed rates) to use in the deterministic cost

PAGE 37

37 node. Next, the process calculations and measurements must be considered. These, in addition to tool life, are the process limiters. The process calculation/measurement node contains the calculations from models or measurements from experiments used to obtain information on process stability, surface location error, and surface roughness. The uncertainty in these calculations depends on the ir inputs, including the uncertain force model coefficients. The unc ertainty in the results is challenging to determine because milling dynamics represent a nonlinear modeling problem. There is also some uncertainty inherent in using the models for process dynamics; however at this point in the research the uncertainty fro m the models has not been quantified. Finally, the cost node gathers the information from all the other nodes and combines it to calculate the expected cost of the process with the selected parameters. A comparison of revenue and cost is completed to obtai n the information for the profit value node. Figure 31. Decision diagram for the milling process. Design Profit Cost Tool life Force model coefficients Revenue a b f t Machining parameter decisions Tool selection decisions Process calc./meas. Tool path

PAGE 38

38 Figure 31 provides the visual representation of all the nodes in the previously described decision diagram. It is seen that information flows from one node to the following nodes so that the required information is made available throughout the process description. These representations enable the remaining work to be completed. Initial Testing Once the milling decision diagram was defined (Figure 31 ), the first step was to test the effectiveness of the model with an example. Therefore, a numerical study was comp leted to test the initial assumptions and the systems sensitivity to the proposed uncertainties. The approach was to choose four optimization examples for a preselected combination of system dynamics workpiece material, Taylor type tool life model, and c utting force model. The four optimization examples are: deterministic (all uncertainties set to zero), tool life uncertainty only, force model coefficient uncertainty only, and both tool life and force model coefficient uncertainties. The tool point FRF w as obtained by modeling a fixedfree 10 mm diameter, 42 mm long cemented carbide cylindrical beam. The properties used to describe the carbide were 550 GPa modulus of elasticity, 14500 kg/m3 density and a structural (solid) damping constant of 0.0015. T hes e parameters were chosen to provide a simple dynamic model for initial testing. Given the model, stability and surface location error behavior was evaluated The mean value of the tool life was obtained using data from 42 tests presented in [ 73]. In this paper, a 10 mm diameter, four tooth tungsten carbide cutter with a 30 degree helix and TiAlN coating was used to machine SKD61 tool steel, corresponding to ISO CMC No. P03.11 alloy steel. The form of the Taylor type tool life model applied here varies slightly from equation 17 because it includes a dependency on the axial depth of cut, b; s ee equation 31, where T is expressed in minutes, v is in m/min, ft is in mm/tooth, and b is in mm. The coefficient and three exponents were determined from l east squares fitting to the test data in [7 3].

PAGE 39

39 = 1 .9549 106 1 6265 0 1024 0 2837 (3 1) In the numerical example, the machining away of a 100 mm cube of tool steel was considered. Manufacturer recommendations for the operating point were used to obtain a basis for comparison. Because manufacturer recommendations were not available for the exact tool workpiece combination used in [73], a similar tool was selected from Sandvik Coromant, which enab led the use of their recommendations. The equivalent tool with approximately the same coating as the tool specified by [73] has the grade designation GC1630 and the closest available geometry was provided by an R216.3410045AC22N tool [74]. Cutting force coefficients were given as 2395 N/mm2 for the tangential direction ( ktc) and 718 N/mm2 for the radial direction ( krc), while the edge coefficients were not provided and were therefore set to zero for this example. The manufacturer recommended cutting condi tions were a surface speed of 88.5 m/min or 2817 rpm for a tool with a 10 mm diameter, using the 1.18 correction factor suggested due to workpiece hardness, and a 0.027 mm/tooth feed rate [74]. = ( 32) = = ( 33 ) = 2 ( + ) + ( 34) = ( + ) ( 35) For the initial study, the machining parameter decisions were discretized into seven values each as shown in T able 31; seven values for each of the four decisions yields a total of 74 or 2401 possible operating conditions. Since the manufacturer only provided suggestions for two out of the four decision parameters (spindle speed and feed per tooth), a suboptimization was performed on the other two in order to obtain the lowest cost. S tability lobes were constructed at

PAGE 40

40 each available radial depth and these were used to select the best possible axial depth. Once the last two parameters (axial and radial depth) were determined, the tool life was calculated using equation 3 1 and the assoc iated cost was determined using equation 1 8, where tool changing time was assumed to be 4 s, the cost per tool was $114 and the operating cost per minute was $1. The machining times for the down milling operation were obtained using the tool path shown in Figure 32 and defined by equations 32 through 35. Given this information, the cost of the operation was determined to be $808.48. In the subsequent 2401 comparisons, the part revenue was set equal to this cost, such that any positive value of profit in dicate d an improvement over the manufacturer recommended operating conditions. Figure 32. Part path used in the machining of the W = 100 mm cube (viewed through the tool center). The path is repeated at multiple levels (into the page); each with the selected axial depth of cut. In step 1 the tool is placed at the correct radial depth. In step 2 the material is removed and in step 3 the tool returns without cutting along the same path. Finally, step 4 repeats step 1. Table 3 1. Discretized machining parameters values used in the initial numerical st udy. a (mm) ft (mm/tooth) b (mm) (rpm) 2.0 0.01 1.0 23000 2.5 0.03 1.5 25667 3.0 0.05 2.0 28333 3.5 0.07 2.5 31000 4.0 0.09 3.0 33667 4.5 0.11 3.5 36333 5.0 0.15 4.0 39000 2 3 same path as 2 (no cutting) 4 W 1 a

PAGE 41

41 Deterministic Optimization In the deterministic example, the maximum profit operating point was selected using the following constraints: the operating point was stable (no chatter), the absolute value of the surface location error or SLE, was less than or equal to 50 m and the average surface roughness was less than 1 m. See Appendices A and B for model descriptions and Matlab code. Among the available choices, the optimal point yielded a profit of $762.23 and was obtained from the parameter set { a, ft, b, } of {4.5 mm, 0.15 mm/tooth, 2 mm, 36333 r pm}. Figure 33 shows the test grid for the optimal radial depth and feed per tooth combination as a function of spindle speed and axial depth. The stability boundary ( the U shaped lobes) and profit contours are also shown (the labels on each contour are e xpressed in dollar values). Filled circles in the grid indicate infeasible points, open diamonds indicate feasible nonoptimal points and the open star identifies the optimum feasible point. Figure 34 shows the SLE contours, rather than the profit contour s where only contours between 50 m are shown Figure 33. Stability diagram for the deterministic example with the lines of constant profit superimposed ( a = 4.5 mm and ft = 0.15 mm/tooth). 554 554 55404704704 7047347347347547547547647647647697697697774774774779779779784784784 (rpm)b (mm) 2 2.5 3 3.5 4 x 104 0 1 2 3 4 5

PAGE 42

42 Figure 34. SLE contours superimposed on the stability diagram for the same case as Figure 3 3. Note that the optimum is in a location with high sensitivity of SLE to small changes in operating parameters (many contours packed together). Uncertain Force Model Optimization For this case tool life was again assumed to be deterministic, but the force model coefficients ( ktc and krc) were considered as uncertain. The HighBaseLow values of the coefficients coefficients ( FHigh, FBase, and FLow, respectively) were obtained by direct encoding assessment of a representative user. For each of the High Base Low force model coefficient values, the optimal profit was calculated. After assuming a symmetric distribution about the base ktc value for aluminium (700 N/mm2 is typical for aluminium alloys), an increment of 150 N/mm2 (a 21.4% change) was selected to correspond to the low and high values. This was then linearly scaled to the 2395 N/mm2 coefficient selected for these examples. The same scaling was applied to the radial coefficient. The approach scaled the coeffi cients for aluminium because the subject of the encoding assessment had more experience cutting aluminium than tool steel. Given the models selected for this case study, uncertainty in the force model coefficients translate d into uncertainty in the results of both the stability and SLE calculations. Higher coefficients tend to reduce the stability limit and increase the absolute value for SLE at any given -36-29-29-21-21-14-14-14-7-7-7-7-7-700000007777777141414141414141414142121212121212121212929292929292929292936363636363636364343434343434343505050505050 (rpm)b (mm) 2 2.5 3 3.5 4 x 104 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5

PAGE 43

43 point. Stability and SLE calculations were repeated three times, once for each encoded value of the fo rce model coefficients. For each of the 7203 (2401 3) combinations, a multiplier value, ( ), was assigned, where F is the degree of the force model coefficient (low, base or high) and S is the setting (the combination of values from the machining para meter decisions). The multiplier value was set to 1 if all constraints were satisfied and 0 if any one of them failed. Then, combining the information from the profit calculations, the selected encoded distributions and the multipliers, an expected profit value was obtained (profit now includes uncertainty so the new result is a distribution). The expected value of each of the 7203 points was calculated using equation 3 6 Note that the profit value used in the equation is the same as the deterministic valu e because the force model coefficients do not affect the cost or revenue calculations directly. Only tool life affects the cost equation directly and, in this case, tool life is deterministic. [ ( )] = ( ) 0 .25 ( ) + 0 .5 ( ) + 0 .25 ( 36) The optimum profit of $751.45 for this case was obtained for parameters {a = 3 mm, 0.15 mm/tooth, b= 2.5 mm, 36333 rpm}. The stability lobe diagram ( a = 3 mm) with variations due to coefficient levels combined with profit cur ves is shown in Figure 35. The upper stability limit corresponds to the low force coefficients and 14 feasible points are shown (diamonds). However, a number of these points are stable for only the base and/or low force coefficient values. For example, the lower left point is unstable for the high force coefficient values. This leads to an M ( Fhigh,S) value of 0 in equation 3 6 and a significantly reduced expected profit at this potential operating point. Note that despite the introduction of new uncertaint ies the optimal decision did not change significantly; this depends on the objective function or value node which is essentially flat over the domain used for this example.

PAGE 44

44 Figure 35. Stability lobe diagram for the uncertain force model optimization at a radial depth of 3 mm and a feed per tooth of 0.15 mm/tooth. Three levels of the stability lobe are shown: the upper dotted line corresponds to the low force model coefficient set, the solid line corresponds to the base coefficient and the high coefficient is the lower dotted line. Uncertain Tool Life Optimization In this case, tool life wa s considered as uncertain and the force model coefficients we re assumed to be deterministic. The deterministic value of tool life obtained from equation 31 is consider ed to be the base value in the encoding exercise performed for this uncertainty. This value was then varied by 15.73% of the base value to obtain the low and high values for the tool life. The percentage is obtained by linearly scaling the users beliefs o n tool life with a base value of 30 minutes and a variation of 4 minutes and 43 seconds. Because the force model coefficients are deterministic and the tool life is not, the multiplier and the profit values in equation 3 6 exchange places and equation 3 7 is obtained. [ ( )] = ( ) 0 .25 ( ) + 0 .5 ( ) + 0 .25 ( 37) 54554 55455470470470473473473475475475476476476476976976976977477477 (rpm)b (mm) 2 2.5 3 3.5 4 x 104 0 1 2 3 4 5

PAGE 45

45 Figure 36. Profit to spindle speed relationship for a = 4.5 mm, ft = 0.15 mm/tooth and b = 2 mm; the upper and lower dashed lines represent th e highest and lowest tool life encoded values, the solid line is the mean tool life and the dotted line near the mean is the expected value of profit as calculated using equation 37. Because the tool life uncertainty does not affect the process limiters (stability, SLE, and surface roughness) the feasible points are the same as those in the deterministic case. For this example, the optimum profit was $761.84 at the same operating conditions as the deterministic case. This is much closer to the deterministic example than the force model coefficient uncertainty, which suggests that force model uncertainties may introduce more uncertainty into the profit value than t ool life uncertainty. However t he similarity between these results and the deterministic results is primarily due to the symmetric distribution chosen for the tool life uncertainty. Figure 36 shows a graphic representation of this selection. The deterministic value of profit shown by the solid line is close to the expected profit value with uncertain tool life indicated by the dotted black line. This agreement is a result of the symmetric variation in profit due to the assumption of symmetric tool life variation. Another tool life st udy [75] suggest s that an asymmetric distribution may better describe actual behavior, which would create a more significant difference between the deterministic and expected values of profit This sensitivity to the low basehigh value selection may be considered in future studies. 0 1 2 3 4 5 x 104 750 755 760 765 770 775 780 785 790 (rpm)P ($)

PAGE 46

46 Uncertain Force Model and Tool Life Optimization This example considered both tool life and force model coefficients to be uncertain. The same encoded values shown in the previous examples were applied. However since neither th e multipliers nor the profit values are deterministic in this case, they cannot be factored out in the profit equation. The generic form is shown in equation 3 8. [ ( )] = 0 .25 ( ) + 0 .5 ( ) + 0 .25 ;klj ; j a 0 .25 ( ) + 0 .5 ( ) + 0 .25 ( 38) The optimum profit in this case ($750.95) was slightly lower than the one calculated in the uncertain force model example, as expected from combining both sources of uncertainty. It was identified at the same parameters as the uncertain force model example. This may not have occurred had the example not been discretized, but it was expected here. Even with both sources of uncerta inty present and a limited domain of 2401 operating conditions, the final result was still a 93% improvement over the m anufacturers recommendations. Recall that the cost of the operation based on the manufacturers recommendation was $808.48 and the part revenue was set equal to this cost, such that any positive values of profit indicate an improvement over the manufacturer rec ommended operating conditions. Table 32 presents all the results of the discrete optimizations for comparison. Table 3 2. Summary of optimization results for comparison. Optimization a (mm) f t (mm/tooth) b (mm) (rpm) Profit ($) Deterministic 4.5 0.15 2.0 36333 $762.23 Tool life uncertain 4.5 0.15 2.0 36333 $761.84 Force model uncertain 3.0 0.15 2.5 36333 $751.45 Tool life and force model uncertain 3.0 0.15 2.5 36333 $750.95 Figure 37 helps explore the optimum profit in the example with both uncertainties present. All parts of the figure show the area in the immediate vicinity of the optimum, the plots are

PAGE 47

47 generated by fixing two of the decisions and plotting profit over the remaining pair. For figures A and B the increase i n axial and radial depths eventually leads to stability or SLE uncertainties penalizing the profit (a similar effect occurs in Figure 37 B below 2 mm due to SLE penalties). However, for the chosen discretization there is no penalty near the optimum point for an increase in feed per tooth. The figure shows a monotonic increase all the way until the maximum selected value of feed per tooth, although the profit increase seems to be reaching a plateau. On each plot the highest point is always the selected optimum of $750.95 and points with $0.00 profit are simply infeasible. A B C Figure 37. Sensitivity to the different parameter choices with respect to spindle speed in the vicinity of the optimum for the uncertain tool life and force model example; A) shows the radial depth, B) shows the axial depth; and C) shows the feed per tooth. 2 2.5 3 3.5 4 x 104 2 3 4 5 x 10-3 0 200 400 600 800 Spindle speed (rpm) Radial depth (m) Profit ($) 2 2.5 3 3.5 4 x 104 1 2 3 4 x 10-3 0 200 400 600 800 Spindle speed (rpm) Axial depth (m) Profit ($) 2 2.5 3 3.5 4 x 104 0 0.5 1 1.5 x 10-4 0 200 400 600 800 Spindle speed (rpm) Feed per tooth (m/tooth) Profit ($)

PAGE 48

48 C HAPTER 4 BAYESIAN UPDATING OF MILLING PROCESS LIMI TERS Improving Knowledge about the Uncertain Milling System Once input uncertainties and their effects have been evaluated, a natural next step is to consider gaining more information about the given system in order to reduce uncertainty and improve the final result. One approach to obtaining new information is to perform experiments and use the data to update knowledge of the system. An important objective for milling research is to reduce uncertainty in st ability predictions; therefore, techniques to express the uncertainty about the stability of system at a given operating condition were explored. Traditionally, the stability boundary is considered to be a step function such that operating points above t he boundary are unstable while those below are stable. However, due to uncertainties in the force model coefficients and to a lesser exten t the FRF, an infinite number of stability boundaries can be calculated, each with a different probability of occurr ence. The steps required to address this situation are to: 1) identify a method that determines the uncertainty of stability information at any given point in the selected domain; and 2) reduce the initial uncertainty using experiments. As shown throughout Chapter 3, the milling optimization problem may be described as a decision problem with uncertainty in not only the force model but also the tool life. The refore, another objective of this research is to calculate the optimal settings and determine the value of perfect information on the resulting uncertainties considering both stability and tool life The value of perfect information provides an upper limit on the usefulness, or value, of any information gathering activity regarding a particular uncerta inty. Note that the value of perfect information is unattainable because new information always includes some level of uncertainty. Uncertainty about the location of the stability boundary is considered first. Once the value of information for the locatio n of the stability boundary has been calculated, the value of an

PAGE 49

49 experiment that determines whether the stability boundary is above or below for a given axial depth and spindle speed combination can be determined. The value of experiment is always lower th an the value of information because it only provides a refinement on the location of the boundary and not the exact location of the boundary. For example, telling a person in the middle of a warehouse that the item he /she is searching for is ahead is less useful than providing the exact location (row number, shelf height, etc.) of the item but it still offers some value I n order to determine the desirability of a stability experiment, the value of the experiment is compared to the value of information to calculate how much value would be obtained by performing the experiment. This information may help answer other questions when used correctly, such as: Where do I test next? and When do I stop testing and start machining? These questions are also addre ssed in this research. Quantifying Stability Uncertainty Quantifying uncertainty in the stability boundary is the first step in determining the associated value of information. However, identifying the milling stability uncertainty is notably different from obtaining the distributions of the uncertainties for the analysis inputs, such as the distribution of force model coefficients. Unlike force model coefficients and the tool holder spindle machine FRF, which are determined from measurements, the stability boundary is obtained from calculations based on a number of inputs (FRF, force model, milling parameters, and tooling information) and the milling model. The output uncertainty for the stability m odel depends directly upon the input uncertainties, their propagation through the model, and the selected operating point. The stability model used here is described in [ 6, 1 4]. It is a frequency domain, analytical stability analysis that allows the direct use of the measured FRF as one of the inputs (no modal fitting). The other inputs for the calculations are the radial depth of cut (used to determine entry

PAGE 50

50 and exit angles for the cutter teeth during material removal) and the number of teeth (both are det erministic decisions), and the force model. For any one set of inputs, the output of this algorithm is the typical stability diagram, which displays the maximum allowable depth of cut without chatter as a function of spindle speed, i.e., the stability boundary. The traditional Taylor series expansion method for uncertainty evaluation [76] is not easily applied in this situation; therefore, a Monte Carlo approach was used to calculate the output distributions. Monte Carlo simulation requires that the force m odel coefficients and FRF uncertainties be described by their mean values and distributions. The methods for describing these distributions vary according to sources of information and the nature of the quantity. In this exercise, force model coefficient u ncertainties were determined through encoding (that same values were used as in the numerical example in Chapter 3). However, due to the FRF s bivariate nature (it is composed of real and imaginary parts, or, alternately, magnitude and phase), obtaining the uncertainty is more comp licated. A previous study [ 72] identified the main contributors of uncertainty in FRF data obtained using an impact hammer and accelerometer (or other linear transducer) as the calibration coefficients of the instruments. Table 41. Gaussian distribution parameters used in the uncertainty comparison between force model coefficients and FRF calibration coefficients. Quantity Mean Standard deviation Units Force model coefficient (radial) 900 100 N/mm 2 Force model coefficient (tangential) 3000 333.3 N/mm 2 Hammer coefficient 926 0.009 N/V Accelerometer coefficient 1027 0.0033 m/s 2 /V As noted, the influence of these uncertainties (FRF and force model coefficients) on the stability diagram was evaluated using Monte Carlo simulation. The input uncertainties were modeled using Gaussian distributions with the parameters shown in Table 41. The results of the

PAGE 51

51 comparison are shown in Figure 41. Part A of the figure shows the results for two standard deviations of variation for the FRF uncertainty, while part B shows the force model uncertainty results, again with two standard deviations of uncertainty. Given these results, the FRF uncertainty was neglected for the remainder of the stability study. The simulation used to obtain these results modeled a 50% radial immersion cut produced by a twotooth cutter. The FRF was constructed using equation 41 with a stiffness of = 1 .6 107 N/m, a damping ratio of = 0 .05 and a natural frequency of 2400 Hz or = 2400 2 rad /s, where is the excitation frequency in rad/s. =1 2 2 2+ 2 ( 41) A B Figure 41. Comparison between (A) FRF and (B) force model coefficient uncertainties as established in [72] and the encoding exercise, respectively. In many engineering problems, assuming Gaussian distributions for the uncertainties is reasonable unless other, more specific, information is available. However in decision anal ysis it is important that distributions accurately reflect the users beliefs/ available information. These input distributions may tend away from Gaussian as they may be skewed or finite, for example. In these situations, beta distributions may provide a better representation of the users beliefs. In addition, due to the nonlinear nature of the stability calculations [ 6, 14] there is no guarantee 0.5 1 1.5 2 x 104 0.5 1 1.5 2 2.5 3 3.5 Spindle speed (rpm)Axial depth (mm) 0.5 1 1.5 2 x 104 0.5 1 1.5 2 2.5 3 3.5 4 Spindle speed (rpm)Axial depth (mm)

PAGE 52

52 that the distribution of the lobes will be Gaussian even if the input uncertainty is Gaussian. Therefore, beta distributions are selected to represent the output uncertainties as well. Beta distributions are defined within the finite domain between and ; the minimum and maximum limits of the uncertain data The shape of the distribution is defined by the parameters r and k as shown in equation 42, where B is the scaled Beta function with parameters r and k Also in equa tion 4 2, x is the value between and at which the user wishes to define the value of the probability density function. Using this scaled beta distribution, it is possible to quantify all the uncertainties required for this exercise. Examp les of a continuous (0, 1, 3, 3) distribution and a histogram of 1000 points from a randomly sampled (2000, 4000, 3, 3) beta distribution are shown in Figure 42. ( | ) = ( ) 1( ) 1, < < (4 2) In this study, the force model coefficients we re assumed to be pe rfectly correlated. As a result, only the tangential value was sampled for the Monte Carlo simulation and the radial value was calculated as 0.3 times the tangential coefficient. The remaining simulation inputs are presented in T able 4 2. A B Figure 42. A) Continuous (0, 1, 3, 3) beta distribution and B) a histogram of a randomly sampled (2000 N/mm2, 4000N/mm2, 3, 3) beta distribution for the tangential direction force model coefficient (FMC). 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 Normalized scaleProbability Density 2000 2500 3000 3500 4000 0 20 40 60 80 100 120 Tangential coefficient (N/mm2)Occurences

PAGE 53

53 Table 42. Inputs for the Monte Carlo simulation Parameter Value Units Stiffness 1.6x10 7 N/m Damping ratio 0.05 Natural frequency 2400 Hz Number of teeth 4 Radial immersion 50% Tangential force model coefficient 30001000 N/mm 2 Radial force model coefficient 900300 N/mm 2 The 1000 randomly sampled coefficients (Figure 42B) were then propagated through the stability model to obtain 1000 stability boundaries over the selected domain (Figure 43). These 1000 boundaries provide a distribution of results at each spindle speed. An example histogram of boundary values at 18000 rpm is provided in Figure 44. For this study, the distributions had the same r and k values at each spindle speed, while the maximum and minimum values of the scaled distribution changed. T his may not be the case in general, however; for example, if FRF uncertainty was considered the bivariate distribution may influence the stability boundary distribution. Fitting the histogram in Figure 44A to a scaled beta distribution using maximum likelihood estimates of the parameters yielded a beta distribution with parameters (2 mm, 3.8 mm, 1.76, 2.56). A normalized version of the distribution (0, 1, 1.76, 2.56) is also shown in Figure 44B. The value of the cumulative distribution obtained from the Monte Carlo simulation at any combination of axial depth and spindle speed indicat es the likelihood that the true stability boundary is below the selected axial depth. In other words, it indicates the probability that a particular cut will be unstable given the information available (i.e., the probability of instability). T he likelihood that a particular parameter combination, s, is stable is obtained by subtracting the probability of instability from one. Because t he beta distribution is defined over a finite domain, the remainder of the probability values with in the domain must be defined separately. In this

PAGE 54

54 case, anything above the calculated distribution was designated zero for unstable, and anything below was stable and was assigned a probability of stability, ( | ) of one. The results are shown in Figure 45. Figure 43. Full set of 1000 stability lobes obtained from the random sampling of the force model coefficients propagated through the Monte Carlo simulation. A B Figure 44. A) Histogram of the distribution of the stability limits at 18000 rpm. B) Normalized beta distribution (0, 1, 1.76, 2.56) for comparison. The output stability distributions are asymmetric, while the input force model uncertainty is symmetric demonstrating the nonlinearity of the stability calculations. 0.5 0.75 1 1.25 1.5 1.75 2 x 104 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 Spindle speed (rpm)Axial depth (mm) 2 2.5 3 3.5 4 4.5 0 20 40 60 80 100 120 Axial depth (mm)Occurences 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 Normalized scaleProbability Density

PAGE 55

55 Figure 45. Probability of stability results using the cumulative distributions functi ons at each spindle speed. White areas indicate stable operation, ( | ) = 1 and black areas indicate unstable operation, ( | ) = 0 Calculating Profit at a Given Setting Once the probability of stability, ( | ) was established, it was possible to calculate the profit at any setting. There are two possible conditions for which profit must be calculated. The first is the unstable case, where it was assumed that the profit was zero dollars. This ignores the possibil ity of fixing or reworking a part which was produced using unstable conditions. The second condition is a stable cut. In t his case, two uncertainty sources were considered: stability (this incorporates force model and FRF uncertainty in the stability resul ts) and tool life (which does not directly appear the stability calculations). In Chapter 3 i t was shown that a value of profit can be calculated at any given milling setting given the value of tool life at that setting, ( ) By combining the encoded t ool life values, the profit calculations, and the probability of stability the expectation of profit, ( ) can then be calculated at a particular milling setting. The uncertainty for tool life wa s obtained from the same encoding exercise as shown in the numerical example in Chapter 3 using equation 31. The machining time was calculated using Spindle speed (rpm)Axial depth (mm) 0.5 1 1.5 2 x 104 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

PAGE 56

56 equation 3 2 through 35 (same tool path). Cost values were determined using equation 18. At each operating point, the encoded values of tool life were used to scale the calculated tool life and obtain the low and high values for the cost calculation. These were then incorporated into the profit calculation in equation 4 3. Finally, the expected profit of milling considering tool life uncertainty with setting s was obtained using equation 4 4. ( ) = ( ) ( 43) ( ) = 0 .25 ( ) + 0 .5 ( ) + 0 .25 ( 44) How ever, tool life uncertainty does not affect stability uncertainty directly, it only affects the expectation of profit in the current formulation of the milling system. Therefore, throughout the next sections the tool life was calculated as a deterministic quantity and the value of profit was obtained using equation 43. Milling Optimization and Experimentation i n t he Presence of Stability Uncertainty Value with N o E xperimentation Stability diagrams are calculated for a particular radial depth of cut. Also, the feed per tooth value does not have a direct impact on the diagram [1 4 ]. Therefore, for this optimization formulation only two of the milling parameter decisions are consider ed (the domain includes only axial depth of cut and spindle speed). Having characterized the stability and tool life uncertainties, as well as the profit, at a given setting, a decision analytic model for the milling optimization was derived A simplified decision diagram and tree [24] which represents the model with no experimental information, is shown in Figure 46. In Figure 4 6, the decision node, Milling setting characterizes the two decisions: spindle speed and axial depth of cut The oval, Mil ling stability, represents stability uncertainty at a given setting. The hexagon is the objective function, Profit. The user desires to make the

PAGE 57

57 Figure 46. Decision diagram and generic decision tree for milling without additional information from experimental data. optimal milling setting decision that maximizes the expected value of profit given the uncertaint y about stability. The arrow in the diagram from Milling setting to Milling stability implies that the probability of stability is influenced by the choice of the spindle speed and axial depth of cut. The arrows from Milling setting and Milling stabili ty into profit indicate that both factors influence the expected profit. A generic decision tree is shown at the right of the figure, where Milling setting decisions are represented by a square, and uncertainty of stability is represented by a circle. Assu ming that the profit is zero when the milling conditions are unstable, the optimal profit is calculated using: = { ( ) ( | ) } ( 45) where is the maximum expected value of the profit that can be obtained given the current state of information about uncertaint ies in the stability boundary and tool life and their effect s on profit. If the user had more or better information about the location of the stability boundary, Value with P erfect I nformation Figure 47. Decision diagram and tree for optimization with perfect information. Milling setting Milling Stability Profit Uncertainty of stability Experiment Setting Milling setting Milling Stability Profit Uncertainty of stability Milling Setting Milling setting Milling Stability Profit Milling setting Milling Stability Profit Uncertainty of stability Experiment Setting Uncertainty of stability Experiment Setting Milling setting Milling Stability Profit Milling setting Milling s tability Profit Uncertainty of stability Milling Setting Uncertainty of stability Milling s etting Milling setting Milling Stability Profit Uncertainty of Stability Milling Setting Milling setting Milling s tability Profit Uncertainty of s tability Milling s etting

PAGE 58

58 he/she might select a different operating point. However, given the current state of information, the value in equation 4 5 is the maximum expected profit. The value of perfect information is the maximum amount the user would be willing to pay to get information about the true location of the stability boundary. For a risk neutral decision maker (one who values uncertain deals at their expe cted value), the value of perfect information is equal to the difference between the value of the deal with perfect information and the value of the deal with no information. A risk seeking user would prefer conditions of greater profit with less regard to the uncertainty, whereas a risk adverse user would be more conservative than the risk neutral user. The value of perfect information is also, by definition, the highest achievable value of any information gathering activity because it gives the best infor mation, the exact location of the stability contour in this case, without any uncertainty in the result. The value of any information obtained from experimental data (information that contains uncertainty) will always be lower. In milling, if the stability boundary was exactly known, the profit at a given setting (either zero or ( )) could be exactly stated. This would enable optimization with perfect information about the location of the stability boundary. Figure 47 shows the decision diagram for this type of situation. Note that information only flows outwards from the milling stability node, which indicates perfect information about stability prior to making any parameter decisions. To calculate the value of the deal with perfect information (the de al is made with the hypothetical clairvoyant to obtain the information) about the stability boundary, the optimal profit that can be obtained for each location of the stability boundary is determined and multiplied by the probability that the boundary abov e the test location, The information is multiplied by the probability of stability to account for the users current belief about the stability

PAGE 59

59 boundary or the belief on what the answer the clairvoyant will give The optimal profit with perfect information is calculated using equation 46, where E denotes the expectation of the value within the brackets. = [ max ( ) | ] ( 46) The difference between and is the value of information if a decision maker is risk neutral. Although it is typically infeasible to obtain perfect information, this number is still important as a comparison to other types of information gathering since it represents the highest possible value. It also serves as a base for comparison for less than perfect information, such as experimental results. Value of E xperimentation Experiments are generally performed to obtain partial information about an uncer tainty or system in question. Although perfect information can never be obtained, partial information may still reduce uncertainty and increase profit, for example In addition, the value of experimentation can provide the user with data to determine the f easibility of proceeding with further experiments. For example, if the value of the experiment is less than the cost of the experiment, then the user may wish to continue working with his/her current state of information rather than Figure 48. Summarized decision diagram and decision tree for a single experimental updating procedure. Milling setting Milling Stability Profit Experiment setting Experiment Stability Uncertainty of Stability Updated Uncertainty of Stability Milling Setting Experiment Setting Uncertainty of Stability Milling Setting Updated Uncertainty of Stability Milling setting Milling stability Profit Experiment setting Experiment stability Uncertainty of stability Updated uncertainty of stability Milling setting Experiment setting Uncertainty of stability Milling setting Updated uncertainty of stability

PAGE 60

60 waste time and money performing the experiment. The general decision diagram and decision tree for a single experime nt are shown in Figure 48. Note that the Experiment setting node and Experiment stability node in the decision diagram may contain information from a number of experiments. Bayesian Predictions of Functions Given the ability to express stability uncertainty and obtain a value for experimentation, it is possible to update stability uncertainties based on experimental results. In many instances, there is not enough information available to generate a stability lobe diagram and only minimal information is available regarding the milling system stability. Therefore, as a first step, a prior with a minimum of information about the relationships between stability, spindle speed, and axial depth of cut was assumed. A methodology for updating using Bayesia n prediction methods in the presence of minimal information about the uncertain stability boundary was then derived. Specifically, it was assumed that the only information available to the user was the maximum allowable depth of cut and a general belief that instability is more likely as axial depth increases. There was no prior information regarding the FRF or force model available at the time of experimentation. Future research can address updating with a more informed prior (see Chapter 6). Bayesian pred iction methods differ fundamentally from curve fitting methods in that they consider a set of predictions as opposed to a single prediction. In other words, while traditional methods try to find a particular function, g(x) to model data, Bayesian methods explore the likelihood that any of a multitude of functions will be equal to the particular function that is used to model the data. This also allows the user to include arbitrary function characteristics such as differentiability class, monotonicity, or p eriodicity, as descriptions of their functional

PAGE 61

61 requirements, which allows multiple information sources to be considered in a mathematically rigorous format. Probability assignment over a functional space is equivalent to applying a multivariate distribut ion over infinite dimensions and then simplifying it so it applies only to the measurable quantities (observables) of interest. The most general set of these is the set of all functionals (any operator that takes in a vector or function as input and output s a scalar) relevant to the phenomenon in question. Typically, this can be expressed in a finite number of dimensions. T his set of all observables is denoted and individual observable is denoted as Oi. There are three formats of the probability distrib utions that are of particular interest: the probability density function (PDF), the marginal probability distribution, and the conditional probability distribution [ 77] The PDF, denoted ( ) states the probability that {O1,O2,On} = { o1,o2,on} for a n indexed set of n observables and any values of oi, 1< i < n. The marginal probability distribution gives the probability that observable J = j given no further information about the rest of the variables. Therefore, i n terms of the PDF it is given by equation 4 7, which is the collapse of the PDF upon all of its dimensions except the one in question. Marginal distributions of multiple observables are defined similarly. = ( ) 1 ( 47) Finally, c onditional distributions 1 | + 1 ( 1, ) denote the probability of {Om+1,Om+2,,On} = {om+1,om+2on} given {O1,O2,,Om} = {o1,o2om}, or, given any partial information on some of the observables, the probability of the rest of the observables. Mathematically, the conditional probability distribution is given by the PDF divided by a marginal distribution such as in equation 48.

PAGE 62

62 1 | + 1 ( 1, ) = ( ) + 1 ,, ( + 1) ( 48) Suppose that the set of observables is chosen as the values of g( x ) for all values of x Also assume that the value of g( x ) for an arbitrary x depends only on the value of g in a small neighborhood of x Therefore, i n terms of the conditional probabilities, for any x1< x2< x3< x4< x5 the distribution of g(x3) would be expressed as equation 4 9. (3) | (2), (4)= (3) | ( 1) (2) ( 4) (5) ( 49) Now g( x ) can instead be represented in terms of its derivative with respect to x ( x ,g ) [ 78] The advantage of this representation is that by equation 49 all ( x ,g ) are independent of each other. If it is assumed that ( x ,g ) are generated from distributions in the exponential family, the PDF can be written in terms of a definite integral. If ( x ,g ) is assumed to be normally distributed for every { x g}, the probability density functional takes the form presented in equation 4 10, where ( ) is the diffusion field, ( ) is the drift field, and ( ) is the creation and killing field. ( )( ( ) ) = ( ) ( ) ( ) 2+ ( ) (4 10) Distributions of this form are known as Brownian distributions and they are ideal for predicting and updating when little is known about the function in question. Assuming any information gained through experiment can be expressed in terms of the function in question ( g) or its derivative ( ) any distribution updated from a Brownian distribution will also be a Brownian distribution. If the experiment measures g directly, the creation and killing field is updated. If an experiment m easures the drift field is updated. The diffusion field may be updated in either situation and can be thought of as a measure of uncertainty at {g,x }.

PAGE 63

63 Updating Milling Stability using Brownian Distributions Updating stability behavior in the spindle spe ed and axial depth domain requires choosing a type of distribution that allows mathematical updating of a highly localized area, due to the peaks present in typical stability lobe diagrams. As a result, a Brownian distribution can be used to describe an i nitial knowledge state. This section describes the procedure for updating stability behavior using a numerical example; Table 43 describes the initial parameter choice for the distribution at the initial knowledge state For this example, knowledge of the axial depth (b) at which the stability boundary occurs is updated as a function of the spindle speed ( In addition, it is assumed that the stability limit is known to occur below an axial depth bmax known beforehand. While t his may be an unrealistic assumption in practice the bound provides an intuitive link between the imposed constraint and c(,b) for the current example. Table 4 3. Distribution parameters. Description Parameter name Value Upper bound on stability limit 1.2 mm Diffusion field ( ) 0.00005 Drift field ( ) 0 Creation and killing field ( ) ( ) = < 0 > 0 N ew knowledge for the updating procedure is obtained through stability testing. In milling stability testing, a machinist selects an axial depth spindle speed combination, performs a test cut, and examines the cut surfaces and/or other metrics (such as the sound of the cut, forces generated during cutting, or tool/workpiece displacements) to determine whether the selected operating parameters produced stable or unstable operation. However, f or this numerical evaluation, the testing wa s simulated by comparing the selected test parameters { } to the predefined/refe rence stability limit, via the stability lobe diagram generated using the parameters listed in Table 44.

PAGE 64

64 Table 4 4. Milling system parameters. Description Value Units Tangential cutting coefficient 2500 N/mm2 Radial cutting coefficient 750 N/mm 2 Number of teeth 4 teeth Tool stiffness 5x106 N/m Damping ratio 0.05 Tool natural frequency 2.4 kHz The knowledge gained from a simple stability test can be expressed as an inequality constraint on b( ). If a test is stable, the stability limit must occur above that point (equation 4 11a). If a test is unstable, the stability limit must occur below that point (equation 4 11b). ( ) < if a test is stable ( 411a) ( ) > if a test is unstable ( 411b) This means that stability testing gives information about the function directly and updating of the stability uncertainties is performed by changing the creation and killing field. Equations 4 11a and 411b are then applied by effectively penalizing all b( ) that do not satisfy the results of the test. Therefore, for any number of stability tests the updated creation and killing fields can be described by equation 4 12. ( ) = < 0 > ( ) < ( ) > 0 ( 412) In order to gauge how performing an experiment affects profit, the decisions that are affected by knowledge of the stability limit should be considered. In this study, the decision situation consisted of selecting a spindle speed axial depth combination {op, bop} at which to mill away a cube of material. Table 4 5 lists the parameters used to calculate profit including the dimensions of the cube. Profit for stable operation, P ( op,bop) was calculated using the approach

PAGE 65

65 provided in Chapter 3 in equations 3 1 to 35 which describe the cutting sequence shown in Figure 32. Profit for unstable operating parameters was assumed zero, although the effects of chatter can sometimes be removed by additional finishing operations. For a fixed knowledge state, the probability of stability for an arbitrary spindle speed axial depth combination can be expressed in terms of a marginal distribution ( equation 413). = ( ) (4 13) Given that the process is described as a stochastic process, the Feynman Kac formula [ 79] can be applied in order to describe the process as a partial differential equation. For the case of stability, a diffusion problem was a reasonable choice given minimal information and the presence of uncertainty. The formulation for the problem is present ed in equation 4 14. Table 4 5. Profit calculation parameters. Description Value Units Width of work piece 500 mm Radial depth 5 mm Tool diameter 10 mm Tool change time 4 sec Number of teeth 4 Machining cost per unit time 1 $/min Cost of tool 114 $ ( ) = ( ) < < > 0 (4 14) with boundary conditions ( 0 ) = 0 ( ) = 0 ( 0 ) = ( ) where and are the limits of the probability density function, originally these values are the maximum allowable axial depth and zero respectively. The solution to this differential equation can be expressed as a series in the form shown in equation 4 15, where the coefficients for the series are provided in equation416.

PAGE 66

66 ( ) = sin = 1 where = (4 15) =2 ( ) sin (4 16) Note that this formulation does not have the diffusion, drift, or creation and killing fields expressed the same way as equation 410. The stability testing performed here does not modify the drift field in any way; therefore, it remains at zero and is not included. The diffusion field is expressed in the k constant of equation 415 and also does not change throughout the tests. The creation and killing field is treated somewhat differently than previously described. It is not expressed in the solution except as a change in the limits of integration ( and ) while proceeding through the different spindle speeds. Essentially, at any spindle speed at which a test is performed, the te st limits of the integration change. If a test is stable, then the lower limit of the integration becomes the axial depth of that test. If a test is unstable, the upper limit changes. Finally, these equations have to be solved both as spindle speed increas es and as spindle speed decreases to obtain the desired symmetric effect for stability updating. If solved in only one direction the new test simply acts as a discontinuity in the solution and the updated PDF would remain unaffected until after it passed t he test spindle speed. Stability behavior near the peaks in the stability lobe diagram does not behave in this manner; it is typically a nearly symmetric peak. To remedy this issue, the equation must be solved from both directions (increasing and decreasin g spindle speeds) and the results from each solution multiplied on a point by point basis. The product is then normalized such that the cumulative probability goes from zero to one. The results of each step are shown in Figure 4 9 for the third cut in the updating sequence, an unstable cut near a 0.5 mm axial depth at just over 10000 rpm. Note the discontinuity prior to encountering a test boundary in both directions (A and B) The updating process is repeated for

PAGE 67

67 each new test so that the effects are cumu lative. The data near 5000 rpm shows the updating results of the two previous test cuts as well, one stable and one unstable cut around a 0.25 mm axial depth. These plots show probability density, and not probability; therefore the data goes from zero to some value and then back to zero with increasing axial depth. Note that because any updating causes a truncation in the PDF and its domain, the resulting distribution has much higher density values than the ones surrounding it (suggested by the significan t color gradient near the update points) such that when changed to a cumulative distribution the probability goes from zero to one from to Figure 49. Probability density functions (each spindle speed has a different PDF) for the numerical solution of the stability updating problem. The three plots represent the different stages of the solution: A) is the solution with increasing spindle speeds, B) is the solution with decreasing spindle speeds, and C) combined and nor malized solution. Using the prescribed algorithm for an arbitrary test history ( equation 4 17), the expected profit of a production run, can be calculated using the probability tree shown in Figure 410. The result is given in equati on 4 18. : = { } 1 < < ( 417) Axial depth (mm) 0.5 0.75 1 1.25 1.5 1.75 2 x 104 0 0.5 1 Axial depth (mm) 0.5 0.75 1 1.25 1.5 1.75 2 2 x 104 0 0.5 1 Axial depth (mm) Spindle speed (rpm) 0.5 0.75 1 1.25 1.5 1.75 2 2 x 104 0 0.5 1 A B C

PAGE 68

68 Figure 410. Probability tree for calculating the expected profit of a production run Pexp op, bop; T ; = ; ( 418) For this example, the operating conditions are chosen to maximize the expected profit, or the value of performing the operation given an arbitrary test history, { } The expected profit is defined using equation 419. { } : = sup, ; (4 19) For a stability test at an arbitrary { t,bt} selection, the probability of each outcome, as well as the expected value of each outcome can be calculated bef ore actually knowing the outcome of the test, using the marginal distribution in equation 413 and the expected value from equations. 418 and 419. Using the probability tree shown in Figure 411, the expected value of the operation after the experiments result may be calculated. Subtracting the value of the profit prior to the result of the experiment being known yields { } the expected value gained from performing that experiment ( equation 420), where is the current knowledge state, denotes augmented with an additional row vector { } and denotes augmented with an additional row vector { } { } = 1 ; { } + ; { } (34) 1 ; ; Value of outcome 0 Outcome Probability of outcome

PAGE 69

69 Figure 411. Probability tree for calculating the value of performing an experiment, { } These probability trees can be chained to represent a series of tests. By using these larger probability trees, the value of an arbitrary number of experiments can be calculated. Therefore, it is possible to assign a value to any arbitrary sequence of experiments, which can then be maximized to determine the optimal testing policy. Since stability tests can be chosen sequentially (i.e., the results of the first test is known when choosing the testing parameters for the second test), the number of test parameters grows exponentially with the number of tests in the sequence (because th e parameters for the second test depend on the result of the first test). If the entire tree is examined, akin to examining all possible occurrences before performing any tests, the testing sequence is planned out ahead of time according to the one that yi elds the highest value. This is called a dynamic search. On the other hand, if the only trees examined are those directly related to the test and the results of a previous test, the user is examining only the part of the tree which maximizes value at each test selection independently. This is called a sequential, or greedy, search. To reduce the computational complexity of this optimization, a greedy strategy was used such that each test, I was chosen by maximizing { } Using this method, a s equence of twelve experiments were chosen and tested against the reference stability lobe diagram. Figure 412 shows the progression of { ; } (the probability that { b} is stable), as the results of a sequence of greedy experiments are sequentially incorporated into the { } Probability of outcome Outcome Value of outcome { ; } 1 { ; } { } { }

PAGE 70

70 prediction. Part A of Figure 4 12 represents the probability distributions prior to any tests being performed showing the same distribution at each spindle speed. The darker areas on the plot indicate areas with a higher probab ility for an unstable cut, the white line on the figure represents the reference stability lobe used to determine the results of each test. Figure 4 12B presents the updated probability distribution after performing four test cuts (two near 5000 rpm, one a t 10000 rpm and one just below 15000 rpm), most of which provided an unstable result and therefore pulled down the black area of the distribution at and near the test locations. After eight tests (Figure 4 12C), several of the tests were stable and they pulled up the white area of the distribution indicating a high Note how the chosen updating functions closely resemble the peaks between the stability functions. Finally, in Figure 412D, after 12 tests, the series nearly match the peak of the stability lobes at around 18000 rpm. Figure 412. Representation of the updating of the probability of stability in a series of 12 tests after A) zero tests, B) 4 tests, C) 8 tests, and D) 12 tests. AAxial depth (mm) 0.5 1 1.5 2 x 104 0 0.2 0.4 0.6 0.8 1 1.2 B 0.5 1 1.5 2 x 104 0 0.2 0.4 0.6 0.8 1 1.2 C Spindle speed (rpm)Axial depth (mm) 0.5 1 1.5 2 x 104 0 0.2 0.4 0.6 0.8 1 D Spindle speed (rpm) 0.5 1 1.5 2 x 104 0 0.2 0.4 0.6 0.8 1

PAGE 71

71 CHAPTER 5 SUPER DIAGRAM Concept A primary challenge in transferring a significant amount of information to a user is providing it in an effective format. Sometimes a user may want to know more about the entire milling parameter space than the point of optimal profit because that optimum may be di fficult or unrealistic to achieve due to factors that were not considered. Similarly, the user may wish to know why a particular area is infeasible. A number of research studies have been completed to observe and address the process limiters that affect th e answers to these questions individually. An empirical tool life model was developed by Talyor in the early 1900s [ 22]. More r ecently, additional insight into tool wear and other machining phenomena w ere detailed by Trent and Wright [ 80 ]. Altintas and Tlu sty also present information regarding the machining process [ 6, 21]. A recent study [ 4] presents a comparison of the different process limiters and their effects on the milling process using a modern milling machine. The results showed that errors due to machining forces and vibrations contributed far more to part errors in highspeed machining than errors caused by temperature changes, CNC contouring error, and quasi static positioning errors. They showed surface location error of over 100 m while the other errors were typically less than 10 m. This is mainly due to the extensive research undergone to minimize errors from sources not dependent on the cutting dynamics (e.g., the machine used in the study incorporated thermal error compensation in its con troller). Therefore, the primary focus of this research is to present users with a clearer picture of errors and limitations caused by the dynamics of the milling system. Each machining setup has a particular machine spindle holder tool part fixture combination; this combination defines the system frequency response function ( FRF ) and force

PAGE 72

72 model. The combination of these two factors provides information about the vibrations of the system during cutting. Analysis of these vibrations provides information about stability and surface location error (SLE) In this work, frequency domain methods a re used to predict both stability and surface location error. Over there past few decades, several frequency domain methods for predicting stability were developed. The basis for these methods was developed by [ 711, 13, 81] and several others. This work uses the method devised in [14], which presents an analytical solution that transforms the time dependent milling equations into a time invariant solution that depends on the selected milling parameters. As noted previously, the data from these analyses is presented in the form of stability lobe diagram s in which stable spindle speed and axial depth pairs appear below the lobes (i.e., the stability limit) and unstable pai rs occur above the lobes. Recently, Schmitz and Mann [19] developed a method to use the same input information to calculate surface location error. Similar to stability, SLE can also be displayed in the spindle speed axial depth domain. However, because SL E is not binary (stable or unstable) the information is typically represented A B Figure 51. A) Stability and B) SLE diagrams (shown in m in the contours) for a sample system. Note that the areas w here high axial depths are achieved (at the lobe peaks) also correspond to the areas of high SLE sensitivity 0.5 1 1.5 2 x 104 0 1 2 3 4 5 6 Spindle speed (rpm) Axial depth (mm) Spindle speed (rpm)Axial Depth (mm) 0.5 1 1.5 2 x 104 1 2 3 4 5 6 100 200 300 400 500

PAGE 73

73 in the form of a contour plot. Figure 511 provides an example of the stability and SLE information presented in separate diagrams for the same sy stem settings and properties. In this work a new user friendly super diagram is introduced which combines the information available from both the stability and SLE predictions. In essence, the well known stability diagram was augmented by adding SLE inf ormation corresponding to a user selected limit. The diagram separates feasible areas where cuts are stable and within the SLE tolerance, from infeasible areas (stable cuts that are out of tolerance and unstable cuts). This research does not define a new prediction method; in fact, the diagram can be developed with any method to obtain the stability and SLE information (temporal finite element analysis for example [1 7, 18, 82]). Rather, the focus is to provide a new approach for viewing all the available i nformation in a single diagram (Figure 5 2) Figure 52. Sample super diagram created using the data in F igure 51 with an SLE user limit of m. Note that SLE blocks out the areas within the stable lobe peaks in this case. Numerical Example For t he super diagram, the predicted stability limit is fully defined for each spindle speed once the system dynamics, radial depth of cut, and force model are specified. However, the Spindle speed (rpm)Axial depth (mm) 0.5 1 1.5 2 x 104 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6

PAGE 74

74 threshold for the acceptable surface location error level is user defined. Th erefore, different users can obtain a different diagram for the same milling system. To introduce the diagram, a typical cutting system was defined and a numerical example was completed. Table 51 provides the cutting system parameters. These include the t ool geometry, system dynamics (represented using symmetric, single degree of freedom modal parameters for simplicity), and the cutting parameters. Table 51. System parameters used in the super diagram numerical examples. Parameter Value Units Stiffness 5x10 6 N/m Damping ratio 0.05 Natural frequency 2400 Hz Tool diameter 6.35 mm Helix angle 45 deg Number of teeth 4 teeth Tangential direction cutting coefficient 700 N/mm 2 Normal direction cutting coefficient 210 N/mm 2 Feed per tooth ( 2 values in the example ) 0.15 mm/tooth 0.075 Radial depth of cut 3.175 mm The limits of surface location error used in this example are based on the absolute value of the calculated SLE. This is not a strict requirement, however. For example, the user may only need to know if a part is overcut (too much material removed) because the part is then tool small and cannot be corrected by a finishing pass or light pass used to remove a small remaining amount of material Therefore, only a constraint on the overcut limit may be desired. If the part is undercut (too little material removed), then finishing cuts can be performed to obtain accurate dimensions (provided the actual dimensions are known). Finally, the user may also choose different limits for over cutting and undercutting. Prior to performing any calculations, users select their operating conditions and the domain of spindle speeds and axial depths for the diagram. This domain is then discretized into a grid of

PAGE 75

75 finite test points. In order to comple te the combined stability and SLE diagram, the stability boundary is first calculated. The test points are then compared to the stability boundary to determine their feasibility (only stable cuts are considered feasible). Next, SLE is calculated for each o f the feasible points. After stability and SLE are determined, penalties are assigned to each test point. Points that are stable and within SLE tolerance levels are not penalized and are labeled with a zero. Points that are stable, but outside of SLE tolerance levels are penalized by one. All unstable points are penalized by two. Given the penalty values, it is possible to generate a corresponding contour plot. Example diagrams are shown in Figure 53 with varying SLE limits (30 m and 70 m for Figures 53A and 53B, respectively) using the information presented in Table 51. This aids in illustrating the effects of user preferences/requirements on the usable machining parameter space. In Figure 53, the black regions are unstable, the grey regions violate the SLE limit, and the white areas are the preferred zones for operation. Note that the grey areas are smaller in Figure 53B due to the relaxed SLE constraint. However, the black area has not changed because stability is not defined by the user, but by the system dynamics. A B Figure 53. Effects of the users SLE limits on the feasible machining domain for errors A) greater than 30 m and B) greater than 70 m at a 0.15 mm/tooth feed rate. Spindle speed (rpm)Axial depth (mm) 1.5 2 2.5 3 3.5 4 4.5 x 104 1 2 3 4 5 6 Spindle speed (rpm)Axial depth (mm) 1.5 2 2.5 3 3.5 4 4.5 x 104 1 2 3 4 5 6

PAGE 76

76 In addition to the user defined SLE limit, the feed rate during the machining process also affects the diagram. Although feed rate has only a second order affect on stability, it is an important parameter in the SLE calculations. For the linear force model used in this case study, the force amplitude scales proportionally with feed rate. In other words, doubling the feed rate doubles the force amplitude. This also doubles the forced vibration amplitude, which doubles the surface location error. Figure 54 displays the same SLE limits as Figure 53, but with half the feed rate. Note that Figure 54B now shows no SLE in excess of the users limit. A B Figure 54. Effects of the users SLE limits on the feasible machining domain for errors A) greater than 30 m and B) greater than 70 m at a 0.075 mm/tooth feed rate. The previous plots show that, in general, SLE varies with spindle speed for a particular axial depth of cut. The source of this SLE behavior is the variation in FRF magnitude and phase with the for cing frequency. Due to the change in phase with tooth passing frequency (spindle speed), the time lag between the force and vibration varies. Therefore, the location of the cutter in its vibration cycle when leaving the surface depends on the selected spin dle speed. The dependence of SLE on the phase lag between the forcing function and displacement causes significant variation near the natural frequency (considering a single degree of freedom system for simplicity) because, for the lowly damped tool point frequency responses typically observed in practice, the phase changes rapidly in this frequency range. Spindle speed (rpm)Axial depth (mm) 1.5 2 2.5 3 3.5 4 4.5 x 104 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 Spindle speed (rpm)Axial depth (mm) 1.5 2 2.5 3 3.5 4 4.5 x 104 1 2 3 4 5 6

PAGE 77

77 The variation in SLE with spindle speed for an axial depth of 4 mm is demonstrated in Figure 55. The significant variation in SLE for spindle speeds is seen near 36000 rpm. This corresponds to a tooth passing frequency of 4 36000 60 = 2400 Hz, which is equal to the selected natural freq uency for the system dynamics (T able 51). This spindle speed is the best speed normally selected for high speed milling applications due to the increased stability limit. Unfortunately, there is a tradeoff between accuracy and the increased material removal rate. L ocal SLE variation is also observed at integer fractions of this best speed. Figure 55. SLE variation with respect to spindle speed at an axial depth of 4 mm and a feed rate of 0.15 mm/tooth. Note the correspondence between sharp changes in SLE with the location of the stable zone peaks in Fig ure 54. Another factor to consider is the effect of the helical cutting teeth geometry on the forces and SLE. The first condition that is considered here is when the combination of system parameters and helical endmill geometry leads to a constant (or nearly constant) cutting force condition. Depending on the users SLE limits, this may expose a feasible zone that would not be present for a straight flute cutter under the same conditions. The location and viability of this parameter space depends on the helix angle and number of teeth on the cutter (the tool geometry 1.5 2 2.5 3 3.5 4 4.5 x 104 -50 -30 -10 10 30 50 50 Spindle speed (rpm)SLE (m)

PAGE 78

78 determines the axial depth at which the phenomenon occurs), as well as the system dynamics (these determine if the axial depth is feasible or not at any given spindle speed). This phenomenon is demonstrated in Figure 56 usi ng the parameters described in T able 5 1. The users SLE limit is 42 m. The axial depth, b, at which constant force occurs is determined using equation 5 1, where d is the tool diameter, is the angular tooth spacing, and is the helix angle of the cut ter. The resulting constant force depth occurs at approximately 4.98 mm, and its integer multiples, using the selected cutter geometry. = 2 tan (5 1) This results in a horizontal step of nearly constant SLE behavior at an axial depth of approximately 4.98 mm (at all stable spindle speeds). For constant force, the deflection calculation simplifies to force divided by the static stiffness (Figure 56). Also, because the force behavior depends on the selected axial depth for helical endmills, SLE varies strongly with axial depth (rather than increasing monotonically as might be expected). This is shown in Figure 5 7 for a fixed spindle speed of 36000 rpm. It is observed that the SLE increases with axial depth up to about 1.25 mm, where it begins to decrease due to the combined effect of phase change and the approach towards the 4.98 mm constant force depth. The slope of sinusoida l variation is caused by the linear scaling of force (and deflection) with axial depth of cut. Figure 58 shows the effect of the interaction between these two phenomenon at the first constant force axial depth. Note that, although the SLE still varies with spindle speed, the variation is much less than in Figure 55. These small variations may be a result of the discretization of the tool and parameter space in the SLE algorithm as well as truncation errors of the axial depth at which these were calculated

PAGE 79

79 Figure 56. Demonstration of the constant force phenomenon. The step behavior is observed at an axial depth of 4.98 mm with a 42 m SLE limit (0.15 mm feed/tooth). Figure 57. SLE variation with respect to axial depth at 36000 rpm and a feed per tooth of 0.15 mm/tooth. Spindle speed (rpm)Axial depth (mm) 1.5 2 2.5 3 3.5 4 4.5 x 104 1 2 3 4 5 6 0 1 2 3 4 5 6 -60 -40 -20 0 20 40 60 80 100 120 Axial depth (mm)SLE (m)

PAGE 80

80 Figure 58. Variation of SLE with respect to spindle speed at an axial depth of 4.98 micrometers and a feed per tooth of 0.15 mm/tooth. Experimental Verification The verification of the super diagram was carried out in two stages, one setup for SLE and another for stability. In both of these setups the part was affixed to a rigid platform (a vice in the SLE experiments and a cutting force dynamometer in the stability experiments) with a compliant tooling system, realized usi ng a long collet type holder and a stubby carbide tool as shown in Figure 59. As a result, the changes in the mass of the part as material was removed did not affect the dynamics that help define the SLE and stability behavior as they would have on a flex ure system with a rigid tool and a flexible part. However, using compliant tooling and a rigid workpiece requires that the spindle dynamics be considered, in general. For the particular machine used in this research (Mikron UCP Vario 600 with a Steptec spi ndle capable of 20 krpm), previous research [83] indicates that the spindle dynamics change significantly with spindle speed (Figure 5 10). As spindle speed increases the natural frequency drops and the dynamic stiffness increases. 1.5 2 2.5 3 3.5 4 4.5 x 104 41.5 42 42.5 Spindle speed (rpm)SLE (m)

PAGE 81

81 Figure 59. Tool and holder combination used in the stability and SLE experiments. On the left is a full view of the tool and holder system while the right is a close up of the end to view the tools insertion into the spindle. Figure 510. Example of changes in the spindle dynamics for the Mikron UCP Vario 600 used in this study as spindle speed increases [83]. Note the decreased amplitude and lower natural frequency of the dominant mode (near 1000 Hz) as speed increases. Surface location error tests were still carried out at higher spindle speeds so that observable error levels could be achieved (SLE tends to decrease with spindle speed and the uncertainty associated with machining and measuring test parts required that the absolute SLE be at the tens

PAGE 82

82 of micrometers level or higher). As will be shown, even in the presence of altered dynamics with increased spindle speed, the general shape of the SLE contours remains identifiable. Stability, however, was verified at lower spindle speeds to reduce the influence of the change in spindle speed dynamics with higher spindle speeds. SLE Verification In order to verify SLE, prior to the actual SLE cutting tests, the system FRF and the force model coefficients were obtained to have a baseline set for comparison between experim ent and prediction. The system FRF was obtained at zero spindle speed just before cutting tests were performed and immediately following machining of the blank test piece. The cutting force coefficients were obtained following the SLE tests using a mechanistic approach [3]. The coefficients were obtained at several different spindle speeds and at a single axial depth of cut of 8 mm (an acceptable depth as suggested by preliminary FRF measurements and typical coefficients). Because the changes in the coefficients were small (the tangential coefficients did not vary by more than 10%), mean values were used in the predictions. A summary of the selected system parameters is shown in T able 5 2. The tool was a Niagara Cutter carbide endmill, model number 85675 wit h 12.7 mm (1/2 inch) removed from the shank end so that when fully inserted in the collet holder only the fluted portion of the tool was visible. Table 52. System and operating parameters for SLE tests. Parameter Value Units Radial depth 3.175 mm Feed per tooth 0.15 mm/tooth Tool radius 6.35 mm Number of teeth 4 teeth Helix angle 30 deg Tangential coefficient 632 N/mm 2 Normal coefficient 122 N/mm 2 Tangential edge coefficient 7779 N/mm Normal edge coefficient 6307 N/mm

PAGE 83

83 The part design used in the SLE experiments is shown in Figure 511 A. The material was 6061 aluminum and the geometry consisted of 24 individual test squares for a single setup. This helped to minimize the effects of alignment errors and changes in individual setups on the resul ts of a set of experiments. The part was machined from a solid block of aluminum just prior to the test cuts and was not removed until the entire sequence of tests was complete. The initial block was approximately 192 mm by 142 mm and 26 mm tall; one side was facemilled to obtain a uniform surface and a 25 mm thickness. The squares were then machined such that they were 20 mm by 20 mm by 15 mm and separated by a distance slightly larger than the diameter of the tool used in these experiments. This allowed for finishing passes (all of which were performed as down milling operations) on the square blocks in order to machine them accurately and leave a uniform reference surface for the measurements of the blocks. Figure 511B shows the profile of a single test block with the reference and test cut surfaces. Also, the sides of the squares were machined such that they were parallel to the X and Y axes of the machine so that when producing the surface any gain mismatch errors between axes in the controller (while small on this machine [1]) were minimized. The individual test blocks were small so that a greater number of tests could be obtained in a single setup. However, this limits the radial depth of cut because larger depths removes more material and leaves behind a smaller measurable surface. A radial immersion of 25% was used in order to leave an adequate measurement surface that excluded the cut entrance/exit transients. SLE was calculated by comparing the distance between the test surface and SLE reference s urface to the commanded radial depth (at 25% radial immersion, 3.175 mm). Also, the size of the 24 blocks was recorded to make sure they were consistent by comparing the distances on the SLE reference surface to the dimension surface. Coordinate measuring machine results showed

PAGE 84

84 that the dimensions of the pre machined blocks (prior to the test cuts) were repeatable with mean dimensions of 20.044 mm (standard deviation of 0.005 mm) in the X direction and 20.038 mm (0.002 mm standard deviation) in the Y direction. A B C Figure 511. A) SLE experimental part as it would appear after test cuts, B) profile of the test squares (not to scale), and C) a picture of the finished part on the machine. Once the test blank was prepared (for a particular test), the tool point FRFs were measured in the X and Y directions; example results are presented in Figure 5 12. The FRFs vary with direction; the Y direction shows more content near the dominant frequency than the X direction, but their magnitudes and natural frequencies are nearly the same. Test surface Dimension surfaces SLE reference surface

PAGE 85

85 Figure 52. FRF measurements results (zero spindle speed) used as input to the surface location error model. Once the simulation input data was gathered SLE predictions were made in order to identify regions that would provide high sensitivity and easily recognizable surface location error features. As a result, test cuts were planned at spindle speeds between 8000 and 12000 rpm (fundamental tooth passing frequencies between 533 and 800 Hz to encompass the dominant mode shown in Figure 5 12) and an axial depth of 8 mm. Test cuts were performed for feed motion in both the X and Y directions. The results of these tests are shown in Figures 5 13 and 514. N ote that in both cases the SLE is smaller than predicted and most of the changes also occur at spindle speeds lower than predicted. This is due to the previously described changes in the spindle dynamics with increasing spindle speed. However, the defining features at the peaks of the SLE are still present (this is marked in the figures with double headed arrows). The observed behavior is supported by the previous FRF measurements on this spindle [83], which exhibited increasing dynamic stiffness and reduce d natural frequency for high spindle speeds. 0 1000 2000 3000 4000 -4 -2 0 2 4 x 10-7 X Real (m/N) 0 1000 2000 3000 4000 -6 -4 -2 0 2 x 10-7 Imaginary (m/N) Frequency (Hz) 0 1000 2000 3000 4000 -4 -2 0 2 4 x 10-7 Y 0 1000 2000 3000 4000 -6 -4 -2 0 2 x 10-7 Frequency (Hz)

PAGE 86

86 Figure 513. Measured and predicted SLE in the X direction (width) of the part as a function of spindle speed. The double headed arrows point at similar features in the contours. Figure 514. Measured and predicted SLE in the Y direction (width) of the part as a function of spindle speed. The double headed arrows point at similar features in the contours. 0.75 0.8 0.85 0.9 0.95 1 1.05 1.1 1.15 1.2 1.25 x 104 20 30 40 50 60 70 80 90 100 110 Spindle speed (rpm)SLE ( m) Predicted Measured 0.75 0.8 0.85 0.9 0.95 1 1.05 1.1 1.15 1.2 1.25 x 104 30 40 50 60 70 80 90 100 110 120 Spindle speed (rpm)SLE ( m) Predicted Measured

PAGE 87

87 Stability Verification While the tool part system was nominally the same (same tool model number a nd part material), small differences in the tool cutting edge geometry and material composition or heat treatment can lead to variation in the force model and small changes in the tool insertion length or clamping force for the collet holder can affect the tool point FRF. Therefore, the FRF and force model measurements were repeated for the stability verification. Force model coefficients were obtained prior to performing the cutting tests and any time a tool was damaged and required replacement. FRF measur ements were repeated each time the tool and holder were removed from the spindle to ensure repeatability and a clear baseline for comparison when performing the stability tests. Table 53 presents the system parameters used in the stability model. The mean tool point FRF from a number of measurements is shown in Figure 515. The mean was adequate for stability prediction and comparison to experiment because the FRFs were consistent between setups. Table 53. System and operating parameters for stability te sts. Parameter Value Units Radial depth 6.35 m m Feed per tooth 0.15 mm/tooth Tool radius 6.35 m m Number of teeth 4 t eeth Helix angle 30 d eg Tangential coefficient 613 N/mm 2 Normal coefficient 149 N/mm 2 Tangential edge coefficient 7019 N/mm Normal edge coefficient 5965 N/mm Similar to the SLE tests, the part was fixed to a rigid platform (a three component force dynamometer) during all the tests. However, the radial depth was increased to reduce the axial depths at which unstable cutting was predicted to occur. The physical l imitation was the 16 mm flute length for the cutting tool. For a 25% radial immersion, the critical stability limit, or the

PAGE 88

88 lowest point on all the lobes, was above 10 mm. By doubling the radial immersion to 50%, the critically stability limit was approxim ately halved. Figure 515. Mean tool point FRF used in the stability calculations. The part used in these experiments is shown in Figure 516 A. It was fixed to the dynamometer, as shown in Figure 5 16 B, using machine screws though the holes at the e nds of the part. The top of the part was faced such that the surface was parallel to the face of the dynamometer and flat. Then, in preparation for the test cuts, a 14 mm deep slot was made with the test tool. Next, the side of the slot that was to be mach ined during the subsequent test cut was finished machined with passes 1 mm deep axially and 0.1 mm deep radially to leave a smooth and accurate surface. Both the slot and the clean up passes were performed at 7250 rpm and a feed per tooth of 0.15 mm/tooth. The clean up pass and the test cuts were all performed as down milling operations. After the surface had been cleaned, the tool was moved to its start point 10 mm to the right of the part at a 6.35 mm radial depth (50% radial immersion) and at the selecte d 0 1000 2000 3000 4000 -4 -2 0 2 4 x 10-7 X Real (m/N) 0 1000 2000 3000 4000 -6 -4 -2 0 2 x 10-7 Imaginary (m/N) Frequency (Hz) 0 1000 2000 3000 4000 -4 -2 0 2 4 x 10-7 Y 0 1000 2000 3000 4000 -6 -4 -2 0 2 x 10-7 Frequency (Hz)

PAGE 89

89 axial depth of cut and spindle speed. Finally, a down milling cut was performed using these parameters. The highlighted surface in Figure 5 16 is the side wall of the test cut. The surface texture of the machined wall was used in combination with the Fourier transform of the force data obtained from the dynamometer in order to determine if chatter (unstable behavior) existed. After the first test cut, the surface was again machined using a cleaning pass and the process was repeated for a new spindle speed axial depth of cut combination. For the given geometry, this allowed for up to 14 test cuts on each block. A B Figure 516. A) A schematic of the test piece with the wall of the test cut surface highlighted and B) a test part with several cut tests alr eady performed. The part is mounted on the force dynamometer. While SLE is a simple comparison of dimensions, identification of instability requires additional analysis. In machining, unstable behavior is described by the existence of self excited vibrations. However, forced vibrations are also present, even if self excited vibration behavior is being exhibited. This complicated the identification of chatter. Forced and self excited vibration can be distinguished in the frequency domain since their frequency content differs: forced vibrations appear at the forcing frequency (the tooth passing frequency) and its harmonics (multiples), while self excited vibrations appear near the dominant system natural frequency (which can be identified from the tool point FRF). An example of a stable (forced vibration) cut

PAGE 90

90 is presented in Figure 5 17 showing both the surface left behind by that cut and the Fourier transform of the steady state part of the force signal (excluding the entrance and exit transients). The cut in this figure was performed at 4500 rpm and a depth of 4.5 mm. The surface on the sidewall displays a repeatable geometry and contains marked nearly vertical cusps. The dominant frequencies in the Fourier transform plot are (from left to right) the once pe r revolution or runout frequency (4500 rpm 1 min/60 s = 75 Hz), the tooth passing frequency (4500 rpm 1 min/60 s 4 teeth/1 rev = 300 Hz), and the first and second harmonics of the tooth passing frequency (600 and 900 Hz). The other, smaller magnitude peaks are either higher harmonics of the tooth passing frequency or the once per revolution frequency. A B Figure 517. Sample data for a stable cut at 4500 rpm and 4.5 mm axial depth. A) Microscope image of the surface of the cut, the inside of the green square indicates the sidewall of the cut. B) The Fourier transform of the force data from 20 to 1500 Hz. A clearly unstable cut, on the other hand, exhibits an irregular and heavily marked sidewall surface and had additional frequency content. The extraneous frequency content may be seen simply as a new peak which is near the systems dominant natural frequency (approximately 700 Hz in this case). However, if this new peak is very near one of the forced vibration frequencies, it can be hard to disce rn or they may even overlap. This makes it hard to quantitatively define and identify chatter and leads to a combined qualitative/quantitative analysis which is subject to 200 400 600 800 1000 1200 1400 2 4 6 8 10 12 x 106 Frequency (Hz)Power (N*m/s)

PAGE 91

91 variation in interpretation from one user to another. Figure 518 shows an example of chatter where the chatter frequency has merged with forced vibration content. Here, the sound of the cut changes (another indication of chatter since added frequency content at a particular frequency changes the tone of the cut), but the surface of the c ut remains fairly repeatable although marks are very heavy and probably outside of user required roughness requirements and SLE tolerances. However, if the frequencies do not merge a new hybrid surface texture is obtained. This is the case for the 4400 r pm 9.5 mm axial depth cut in Figure 519, where chatter occurs near 740 Hz and the forced vibration content appears at 730 Hz. Two more examples of clearly defined chatter are presented in Figures 5 20 and 521. An overview of the stability results is pre sented in Figure 5 22. In the figure the red Xs symbolize the unstable cuts previously presented and the green squares are cuts with no evidence of chatter. The black circles, however, indicate cuts that presented some evidence of chatter (nonuniform sur face texture or small peaks at the chatter frequencies identified from deeper A B Figure 518. Sample data for an unstable cut at 4500 rpm and 9.5 mm axial depth of cut where the chatter frequency matches some of the content from the stable cut. A) Microscope image of the surface of the cut, note the crosshatching along the surface. B) The Fourier transform of the force data from 20 to 1500 Hz, the peak at 750Hz would normally correspond to the ninth harmonic of the runout frequency, but it is too large indicating that a chatter frequency expressed itself at or very near this frequency. 200 400 600 800 1000 1200 1400 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 x 107 Frequency (Hz)Power (N*m/s) Chatter

PAGE 92

92 A B Figure 519. Sample data for an unstable cut at 4400 rpm and 9.5 mm axial depth of cut where the chatter frequency is slightly off from the stable content. A) Microscope image of the surface of the cut, note the crosshatching and diagonal marks along the surface. B) The Fourier transform of the force data from 20 to 1500 Hz, the peak at 740Hz is chatter and t he smaller peak at 730 Hz. A B Figure 520. Sample data for an unstable cut at 6250 rpm and 9.5 mm axial depth of cut. A) Microscope image of the surface of the cut, note the uneven crosshatching along the surface. B) The Fourier transform of the force data from 20 to 1500 Hz, the peak at 730Hz is expected and the smaller peak is due to the chatter frequency at 715 Hz. Here the chatter amplified a peak that was supposed to be small. cuts). Note the significant increase in allowable axial depth of cut at 5500 rpm relative to the other test speeds. This demonstrates the usefulness of stability prediction in selected machining conditions at the process planning stage. 200 400 600 800 1000 1200 1400 2 4 6 8 10 12 14 x 106 Frequency (Hz)Power (N*m/s) 200 400 600 800 1000 1200 1400 2 4 6 8 10 12 14 x 106 Frequency (Hz)Power (N*m/s) Chatter Expected content Expected content Chatter

PAGE 93

93 A B Figure 521. Sample data for an unstable cut at 7000 rpm and 9.5 mm axial depth of cut. A) Microscope image of the surface of the cut, note the crosshatching and diagonal marks along the surface. B) The Fourier transform of the force data from 20 to 1500 Hz, the peak at 700Hz is expected and the larger peak at 720 Hz is chatter. Here the chatter simply overpowers the expected frequency data similar to the test at 4400 rpm. Figure 522. Comparison of the results of the stability experiments with the theoretical stabi lity lobes. The blue line represents the modeled stability lobes using the mean force model coefficients and the mean FRF. The crosses, squares, and circles are the unstable, questionably stable/unstable, and stable cuts respectively. 200 400 600 800 1000 1200 1400 1 2 3 4 5 6 7 8 9 10 11 x 106 Frequency (Hz)Power (N*m/s) 4000 4500 5000 5500 6000 6500 7000 7500 0 2 4 6 8 10 12 Spindle speed (rpm)Axial depth (mm) Expected content Chatter

PAGE 94

94 Closer observation of the roughness of two of the test cut conditions reveals additional information relevant to the discrepancies observed in the stability behavior. Sections of the cuts performed at 4400 rpm, and at axial depths of 4.5 mm and 9.5 mm were measured using a scanning white light interferometer. The results are shown in Figures 523 and 524. During stable cutting each tooth engages the surface equally, assuming there is negligible run out between the teeth. So, given that the feed per tooth for these experimental cutting operations was set to 0.15 mm/tooth, feed marks would be expected every 0.15 mm However, the markings present in Figure 523 occur every 0.6 mm indicating that a period of once per revolution or once every 4 tooth engagements This suggests signi ficant run out of the cutter teeth (i.e., different radii). Similarly, the larger markings on the unstable occur at near 0.6 mm at their widest points in Figure 5 24. The direction of the profile in Figures 523B and 524B is consistent with the feed direc tion of both cuts. The images in part A of the figures show the full axial depth of each cut perpendicular to the profile directions and are stitched together from multiple images using the MetroPro software provided with the Z ygo Newview 7200 scanning whi te light interferometer. Although the peak to peak values of the profiles are unexpectedly similar (due to the runout), they are less consistent during unstable cutting, as expected due to the self excited vibrations caused during a chatter event. After experimentation, the runout of the tool was measured a long the positive and negative X and Y axes of the machine; the measurements were consistently near 20 m peak to peak between the four teeth on the cutter as it was rotated in the spindle. The subse quent variation in chip thickness was not considered in the analytical approaches prescribed here and contributes to the small differences between experimental and predicted results.

PAGE 95

95 A B Figure 523. Scanning white light interferometer image of the s urface left during stable cutting at 4400 rpm and 4.5 mm axial depth. Image A) is the color contoured image of stable cutting, with the line indicating the location at which the profile, B), was obtained. Note the peaks occur every 0.6 mm indicating signif icant run out..

PAGE 96

96 A B Figure 524. Scanning white light interferometer image of the surface left during unstable cutting at 4400 rpm and 9.5 mm axial depth. Image A) is the color contoured image of the unstable cut, with the line indicating the location at which the profile, B), was obtained along the feed direction. Note the peaks occur every 0.6 mm indicating significant run out. The incomplete squares occur due to areas outside of the focus range of the machine as this is a stitched image.

PAGE 97

97 C HAPTER 6 CONCLUSIONS Milling Optimization using Decision Analysis Uncertainty in the milling process is present in each of its stages, including process planning (where the tool path is defined based on the desired part geometry), the milling process itself due to the influence of the dynamics, for example, and validation of the final part characteristics (such as dimensions and fatigue characteristics). It is important that the user consider these uncertainties when planning a milling process (i.e., making decisions) because they may affect his/her choice of operating conditions and maybe even their constraints. This research presented a structure for organizing and evaluating milling process uncertainties using Decision Analysis and Bayesian methodologies. The first step in this evaluation was analyzing the milling proc ess to identify the decisions controlled by the users and designers (milling operating parameters, part design and system selection), the process limiters (stability, surface location error, surface roughness, and tool life), the input uncertainties (the m illing system frequency response function, FRF, and the force model for process simulation), and finally the value inherent to the process given all the previous information (profit). After identifying the information categories, the influences of the indi vidual uncertainties were explored and the effects of the dominant force model uncertainty on the stability model were investigated. As a result, a method to quantify the uncertainty in the stability model as a probabilistic quantity was determined. Exper imentation is a simple and rapid way of obtaining partial information about an uncertain system. One typical experiment on the modern shop floor is using test cuts to identify stable milling parameters. In this research, a methodology was implemented to update stability uncertainties using only simple stability information (stable or unstable only). This method

PAGE 98

98 provided an approach for optimal selection of the next test parameter set (spindle speed and axial depth of cut) given the results of previous tests using Bayesian updating methodologies and very limited knowledge of stability behavior. Prior to the stability tests, the only knowledge of stability the user typically has available is that generally the probability of unstable behavior (chatter) incre ases as axial depth increases and that the maximum depth is constrained by the flute length on the cutter. For the study completed here, this general case was followed such that there was no prior knowledge of stability inputs (FRF and force model) to predict the spindle speed dependence of the stability lobe behavior. A numerical example was performed using this methodology to demonstrate its efficacy. Additional Research Possibilities Involving Decision Analysis in Milling Optimization : Similar to stability uncertainty due to imperfect knowledge of the FRF and force model, surface location error (SLE) predictions will also be uncertain. Because the part errors associated with SLE can be significant, it is naturally of interest to quantify and update SLE uncertainty. Also, a useful improvement to the stability updating scheme would be to improve testing value by starting the updating from a more informed/complete prior and including additional sensor information from experimentation. An example of a more co mplete prior would be if the user obtained FRF and force model data to generate a distribution of the stability boundaries and other process limiters and then used that instead of the more general, less informed prior applied in the numerical example. Howe ver, updating using the informed prior might be slightly different given the structure of the informed distributions. Therefore, additional research is required in this area. The end result of all the uncertainty propagations should be: 1) a single domain or a set of domains with continuous probability functions that reflect the users beliefs concerning the

PAGE 99

99 feasibility of all available cutting parameters based on process limiters (tool life, stability, SLE, surface roughness, etc.); and 2) their propagati on through the profit calculation such that a continuous optimization routine can be used over the entire domain to find optimal cutting conditions for the process. Once all the information is gathered and processed, the optimization can then be performed in real time during machining processes, assuming appropriate sensors are available and incorporated into a smart machining algorithm. The viability of a sensor can also be determined using decision analysis by comparing the information obtained using the sensor to the effect this information would have on the expectation of profit in a particular setup. The choice of sensor would then depend on which one would provide the most information at the least cost and effort for a particular setup and budget. Super D iagram Using all the different process limiters at once for the decision analytic portion of this work posed an interesting question: How can all the information from process limiters and uncertainties be presented to the user in a n easy to understand diagram? Based on previous work, the user can base parameter selections on a stability diagram and contour plot of the surface location error over the same spindle speed axial depth of cut parameter space. However, this may be confusing and requires anal ysis and experience, which may be unavailable. In particular, the surface location error plot contains a wealth of information, but much of it is not immediately identifiable from a cursory view. The solution was to identify the SLE points that cause the p art to be out of tolerance and label them differently than those that were within limits, forming a binary information set similar to the traditional stability boundary. These two binary information sets were then combined using a penalty system (feasible points received no penalty, out of tolerance SLE received a penalty of one, and unstable behavior received a penalty of two) to generate a single diagram that

PAGE 100

100 showed both process limiters over the same domain in a grey scale format. This was referred to as the milling super diagram. This diagram is unique because part of it is user defined ; the SLE limit is based on the user input. This means that two users with different tolerance requirements, but otherwise identical systems, would have different diagra ms. The predictive algorithms for both SLE and stability were then verified experimentally to validate the super diagram information. Future Work Involving the Super Diagram : The super diagram is a versatile tool as it stands. However, the diagram can be e xpanded further. In addition to stability and SLE, other process limiters can be added. Tool life and surface roughness predictions are prime candidates for inclusion. Surface roughness for stable cutting depends strongly on the feed per tooth, but runout of the individual cutter teeth (clearly an uncertain variable) also plays an important role. The SLE prediction algorithm could be adapted to determine surface roughness. Tool life, on the other hand, is typically identified empirically and fit using a pow er law. Alternately, the cutting force can be monitored as the tool wears and an empirical model defined to update the cutting force coefficients based on the amount of material removed, for example. These models can be implemented and the results added to the diagram. In addition, uncertainties can also be added to the diagram such that the borders between regions become continuous as opposed to stepped. Finally, an optimization using the super diagram could be performed. Previously, it was noted that ther e are uncertainties in all of predictions. However, quantifying these uncertainties and experimenting to reduce the uncertainties so that the user can obtain an acceptable prediction takes time and money, which reduces profits. An alternative to directly q uantifying these uncertainties is to select safety margin further from the infeasible points and limit the area available for optimization. For example, any point within 100 rpm of the predicted infeasible

PAGE 101

101 boundaries could be excluded from the optimizati on and points a pre selected distance from an infeasible axial depth could also be excluded. The remainder of the feasible domain can then be used in a global optimization scheme such that a more conservative optimal solution can be identified.

PAGE 102

102 APPENDIX A STABILITY AND SLE MODEL DESCRIPTIONS AND CODE. Stability The Altintas and Budak [ 14] milling stability algorithm is used throughout this research. This algorithm expands the time varying coefficients of the milling dynamics equations (which relate milling forces and system vibrations) which depend on the angular orientation of the cutter as it rotates through the cut, into a Fourier series and then truncated to include only the mean term (a constant). The problem is then expressed as an eigenvalue problem The general form selected is = 0 where represents the eigenvalues for the system I represents the identity matrix, and A is an expression summarizing the milling dynamics relevant to stability. The formulation of the A mat rix is shown in Equation A1, where Gx/y and Hx/y are the real and imaginary parts of the tool point FRFs in the x and y directions, respectively. The terms xx, xy, yx, and yy depend on the selected radial immersion the cutting direction (up milling or down milling), and the cutting force coefficient, = as shown in Equations A2A5 = ( + ) + ( + ) + (A1) =1 2 [ cos 2 2 + sin 2 ] (A2) =1 2 [ sin 2 2 + cos 2 ] (A3) =1 2 [ sin 2 + 2 + cos 2 ] (A4) =1 2 [ cos 2 2 sin 2 ] (A5) The resulting eigenvalues are then used to calculate the limiting axial depth of cut between chatter and stable cutting (Equation A6). The spindle speeds that correspond to the axial depths calculated in Equation A6 are obtained using Equation A7. In these equations, blim is the limiting

PAGE 103

103 axial depth of cut without chatter (the stability boundary) c is the chatter frequency ( if chatter occurs ), and j and N is are integer s ( j = 0 1, 2, and N = 1, 2, 3...) corresponding to the lobe number (they are numbered from right to left) and the number of teeth on the cutter The reader may note that these equations differ slightly from those in [11] because the eigenvalue problem formulation has been slightly modifie d. =2 ( ) ( ( )2+ ( )2) (A6 ) = ( + 2 ) where = 2 tan 1 ( ) ( ) (A7 ) The previous equations are used for both eigenvalues (these are typically a complex conjugate pair); the lowest value at each spindle speed defines the actual stability limit. In the algorithm used for this research, the resulting stability lobe s are trimmed to identify the lowest boundary (the theoretical limit between stability and instability) at any give spindle speed since the calculated lobes (from the two eigenvectors) can cross one another. The stability code presented in this appendix contains an additional section which checks the selected test points used in the numerical examples against the calculated stability lobe, to determine their stability condition.

PAGE 104

104 % Ral E. Zapata % Altintas Stability Lobes Trimmed with search algorithm for stable and % unstable individual points. % October 2006 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [Combination_mat,SS_final,Blim_final] = stability_alt_trim(kt, kn, N, phis, phie, spin_speed, Num_lobes, ADOC_vec, SS_vec, f, Gxx, Gyy) % This function provides the full set of tested points along with the % spindle speed and Blim vectors for the stability lobe diagram they were % compared to. The last column on Combination_mat is the stability % multiplier where 1 means that the parameter set is stable and zero means % the parameter set is unstable. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Calculations %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% w_temp = f*2*pi; w = zeros(1,length(w_temp)); for cnt = 1:length(w_temp) w(1,cnt) = w_temp(cnt); end % Enter in the Cutting Coefficient Information Ktc = kt; Krc = kn; % calculate cutting coefficients for monte carlo run Kta=Ktc; Kra=Krc/Ktc; % -----------------------------------------Altintas Code alphaxx=0.5*((cos(2*phie)2*Kra*phie+Kra*sin(2*phie))(cos(2*phis)2*Kra*phis+Kra*sin(2*phis))); alphaxy=0.5*((sin(2*phie)2*phie+Kra*cos(2*phie))( sin(2*phis)2*phis+Kra*cos(2*phis))); alphayx=0.5*((sin(2*phie)+2*phie+Kra*cos(2*phie))( sin(2*phis)+2*phis+Kra*cos(2*phis))); alphayy=0.5*((cos(2*phie)2*Kra*phieKra*sin(2*phie))( cos(2*phis)2*Kra*phisKra*sin(2*phis))); % Calculate and Sort Eigenvalues for cnt=1:length(w) A=[alphaxx*Gxx(cnt) alphaxy*Gyy(cnt);alphayx*Gxx(cnt) alphayy*Gyy(cnt)]; E=eig(A); temp=E(1); eigen1(cnt)=temp; temp=E(2); eigen2(cnt)=temp; if (cnt>1)

PAGE 105

105 dot_prod1=real(eigen2(cnt))*real(eigen2(cnt1))+imag(eigen2(cnt))*imag(eigen2(cnt1)); dot_prod2=real(eigen2(cnt))*real(eigen1(cnt1))+imag(eigen2(cnt))*imag(eigen1(cnt1)); if(dot_prod2>dot_prod1) temp=eigen2(cnt); eigen2(cnt)=eigen1(cnt); eigen1(cnt)=temp; end end end eigen1=eigen1'; eigen2=eigen2'; % Calculate alim values for each eigenvalue alim1=(2*pi/N/Kta)./((real(eigen1)).^2+(imag(eigen1)).^2).*(real(eigen1).*(1+ (imag(eigen1)./real(eigen1)).^2)); alim2=(2*pi/N/Kta)./((real(eigen2)).^2+(imag(eigen2)).^2).*(real(eigen2).*(1+ (imag(eigen2)./real(eigen2)).^2)); % Choose positive Values of alim1 [index1]=find(alim1>0); alim1=alim1(index1); alim1=alim1*1000; w1=w(index1).'; psi1=atan2(imag(eigen1),real(eigen1)); psi1=psi1(index1); epsilon1=(pi2*psi1); % Note that this is a column vector %Choose positive Values of alim2 [index2]=find(alim2>0); alim2=alim2(index2); alim2=alim2*1000; w2=w(index2).'; psi2=atan2(imag(eigen2),real(eigen2)); psi2=psi2(index2); epsilon2=pi2*psi2; % Note that this is a column vector %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Matrix Formulation N1_matrix = zeros(length(epsilon1),Num_lobes); N2_matrix = zeros(length(epsilon2),Num_lobes); for cnt1 = 1:Num_lobes N1_matrix(:,cnt1) = (60/N)*w1./(epsilon1+2*(cnt11)*pi); N2_matrix(:,cnt1) = (60/N)*w2./(epsilon2+2*(cnt11)*pi); end for cnt2=1:length(spin_speed)

PAGE 106

106 SS_final(cnt2)=spin_speed(cnt2); % Initiallize Search Matrices blim1_temp = zeros(length(epsilon1),Num_lobes); blim1_temp(1,:)=1e20; blim2_temp = zeros(length(epsilon2),Num_lobes); blim2_temp(1,:)=1e20; % Populate search matrices by linearly interpolating the desired SS points. for cnt3 = 1:Num_lobes for cnt4=2:length(epsilon1) if(spin_speed(cnt2)>N1_matrix(cnt4,cnt3) & spin_speed(cnt2)N1_matrix(cnt41,cnt3)) blim1_temp(cnt4,cnt3)=alim1(cnt41)+(alim1(cnt4)alim1(cnt41))*(spin_speed(cnt2)N1_matrix(cnt41,cnt3))/(N1_matrix(cnt4,cnt3)N1_matrix(cnt41,cnt3)); else blim1_temp(cnt4,cnt3)=1e20; end end for cnt4=2:length(epsilon2) if(spin_speed(cnt2)>N2_matrix(cnt4,cnt3) & spin_speed(cnt2)N2_matrix(cnt41,cnt3)) blim2_temp(cnt4,cnt3)=alim2(cnt41)+(alim2(cnt4)alim2(cnt41))*(spin_speed(cnt2)N2_matrix(cnt41,cnt3))/(N2_matrix(cnt4,cnt3)N2_matrix(cnt41,cnt3)); else blim2_temp(cnt4,cnt3)=1e20; end end end Blim1_min = min(min(blim1_temp)); % Find minimum boundary for alim1 lobes at this SS Blim2_min = min(min(blim2_temp)); % Find minimum boundary for alim2 lobes at this SS Blim_final(cnt2) = min([Blim1_min Blim2_min]); % Find global minimum boundary for this SS end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Test parameters stability check Counter = 0; Combination_mat = zeros(length(ADOC_vec)*length(SS_vec),3);

PAGE 107

107 for cnt50 = 1:length(ADOC_vec) for cnt51 = 1:length(SS_vec) Counter = Counter+1; index_low = max(find(SS_final < SS_vec(cnt51))); index_high= min(find(SS_final >= SS_vec(cnt51))); SS_left = SS_final(index_low); SS_right = SS_final(index_high); B_lim_left = Blim_final(index_low); B_lim_right = Blim_final(index_high); Blim_value = B_lim_right+(B_lim_leftB_lim_right)*(SS_vec(cnt51)SS_right)/(SS_leftSS_right); if Blim_value > ADOC_vec(cnt50) STAB_mult = 1; elseif Blim_value <= ADOC_vec(cnt50); STAB_mult = 0; end Combination_mat(Counter,:) = [ADOC_vec(cnt50) SS_vec(cnt51) STAB_mult]; end end Combination_mat;

PAGE 108

108 Surface Location E rror This research uses the analytical, frequency domain solution formulation to surface location error developed in [ 19]. There are two assumptions to consider where using this frequency domain formulation. First, the influence of regeneration is neglected b ecause the machining process is stable and governed by steady state vibrations ; therefore, for any given SLE calculation to be valid the cut must be stable. Second, although tool vibrations occur in both the x and y directions, y direction vibrations dominate the surface location ( the selected reference coordinate frame uses the x direction as the feed direction ) since this is the direction in which the error is measured Given these conditions, the solution method can be divided into three steps: 1. Expr ess the y direction cutting force in the frequency domain. There are several methods to approximate this using the force model, but for this research a simulation of a single rotation of the tool (without regeneration considerations, i.e. using an infinite ly rigid tool) is modeled to obtain the force vector as a function of angular location of the tool. Only one rotation is necessary because this can be concatenated as many times as necessary since subsequent rotations will be identical. 2. Transform the force vector to the frequency domain using a fast Fourier transform and determine the y direction displacement in the frequency domain by multiplying the y direction cutting force by the y direction FRF or ( ) ( ) as shown in Equation A8. ( ) = ( ) ( ) ( ) (A 8) 3. Inverse Fourier transform the resulting displacement data and sample at the cut entrance (for up milling) or exit (down milling) to determine surface location error. The SLE code presented here calculates the values on a point by point basis, so it must be executed for each test point. However, the code can be modified such that it groups all spindle speeds at a particular axial depth into one b atch. This is possible because during a particular rotation (if regeneration is neglected) the force vector is independent of spindle speed for a given axial depth if it is expressed as a function of angle and not time. The angular formulation is then

PAGE 109

109 tran sformed for each spindle speed to calculate the appropriate time vector. This is less computationally expensive than calculating the force vector repeatedly.

PAGE 110

110 % T. Schmitz (10/20/05) % This is a program to find the SLE in helical peripheral end milling using a frequency domain approach. % Regeneration is not considered. function [e] = sle_navy_DCcorr_axialselect(fmeas, FRFmeas, kt, kn, kte, kne, m, beta, d, phistart, phiexit, omega, b, ft, ADM) % Consider only up and downmilling cases if phistart == 0 flag = 0 + 2*ADM*tan(beta)/d; % upmilling else flag = 180 + 2*ADM*tan(beta)/d; % downmilling end DTR = pi/180; % degrees to radians conversion % Simulation specifications n = 11; steps = 2^n; % steps for one cutter revolution, int dt = 60/(steps*omega); % integration time step, s dphi = 360/steps; % angular steps size between time steps, deg if beta == 0 % straight teeth db = b; % discretized axial depth, m else % nonzero helix angle db = d*(dphi*DTR)/2/tan(beta*DTR); end steps_axial = round(b/db); % number of steps along tool axis tooth_angle = 0:360/m:(360360/m); % angles of m cutter teeth starting from zero, deg % Initialize vectors teeth = round(tooth_angle/dphi) + 1; phi = linspace(0, (steps1)*dphi, steps); Force_y = zeros(1, steps); %************************** MAIN PROGRAM ****************************** for cnt1 = 1:steps % time steps, s for cnt2 = 1:m teeth(cnt2) = teeth(cnt2) + 1; % index teeth pointer one position (rotate cutter by dphi) if teeth(cnt2) > steps teeth(cnt2) = 1; end end Fy = 0; for cnt3 = 1:m % sum forces over all teeth, N for cnt4 = 1:steps_axial % sum forces along axial depth of helical endmill, N phi_counter = teeth(cnt3) (cnt41); if phi_counter < 1 % helix has wrapped through phi = 0 deg phi_counter = phi_counter + steps; end phia = phi(phi_counter); % angle for given axial disk, deg

PAGE 111

111 if (phia >= phistart) & (phia <= phiexit) % verify that tooth angle is in specified range for current disk, deg h = ft*sin(phia*DTR); % chip thickness, m ftan = kt*db*h + kte*db; frad = kn*db*h + kne*db; else % tooth angle is outside range bounded by radial immersion ftan = 0; frad = 0; end Fy = Fy frad*cos(phia*DTR) + ftan*sin(phia*DTR); % N end % cnt4 loop end % cnt3 loop Force_y(cnt1) = Fy; end % cnt1 loop %************************** END OF MAIN PROGRAM ****************************** % Compute FFT of Fourierbased ydirection force [FY, freq] = spec(Force_y', 1/dt); % compute FFT (in spec.m, there should be no multiplication by T) FY = FY/(2^n); % correct magnitude to N % Define ydirection FRF on proper frequency vector, freq index = find(freq >= min(fmeas) & freq <= max(fmeas)); freqtemp = freq(index); FRFytemp2 = interp1(fmeas, FRFmeas, freqtemp, 'spline'); FRFytemp1 = ones(1, (index(1)1))*FRFytemp2(1); FRFytemp3 = ones(1, (length(freq)index(length(index))))*FRFytemp2(length(FRFytemp2)); FRFy = [FRFytemp1 FRFytemp2' FRFytemp3]'; Yf = FY.*FRFy; % F X/F = X Yf_dc = Yf(1); % DC component Yf(1) = 0; % DC extraction y = real(ifft(Yf*(2^n))); % convert to timedomain y = y+Yf_dc; % DC insertion to create full signal t = 0:2*dt:(length(y)1)*2*dt; % new time vector after inverse FFT % Use automatic method to sample y tfirst = (flag*DTR)/(omega/60*2*pi); % time for first SLE point, s index = find(t > tfirst); first_point = index(1) 1; % first point for SLE index = first_point:round(60/(omega*m)/(2*dt)):length(t); y_sampled = y(index); % sampled position vector -SLE values e = y_sampled(length(y_sampled)); % record SLE, m

PAGE 112

112 APPENDIX B DETERMINISTIC EXAMPLE CODES FOR PROFIT M AXIMIZATION Main Code % CALL ALL with FMC encoding % Raul Zapata % December 2006 clc % close all clear all tic %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % FRF % %%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Load modeled FRF from beam theory % % load rigidfreetool_FRF_data50_1e.mat % % % % fmeas = f; % % FRFy = H11; % % % % FRFy(1) = FRFy(2); % % FRFx = FRFy; load rigidfreetool_FRF_data_42casestudy_1e.mat w = f.'*2*pi; fmeas = f; Gxx = H11; Gyy = H11; Gxx(1) = Gxx(2); Gyy(1) = Gyy(2); FRFx = Gxx; FRFy = Gyy; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % MODAL FIT % % % % % Tool description % % % % ky = [5e6 7e6]; % N/m % % % % zetay = [0.05 0.05]; % % % % wny = [800 900]*2*pi; % rad/s % % % % my = ky./(wny.^2); % kg % % % % cy = 2*zetay.*(my.*ky).^0.5; % Ns/m % % % % % % % % % Define ydirection FRF (could read in measurement from TXF)

PAGE 113

113 % % % % fmeas = 0:0.1:5000; % Hz % % % % w = fmeas*2*pi; % rad/s % % % % FRFy = (wny(1)^2/ky(1))./(wny(1)^2 w.^2 + i*2*zetay(1)*wny(1).*w); % % % % for cnt = 2:length(ky) % % % % FRFy = FRFy + (wny(cnt)^2/ky(cnt))./(wny(cnt)^2 w.^2 + i*2*zetay(cnt)*wny(cnt).*w); % % % % end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % .TXF FILE CALL % % % % % Assumptions: % % % % % 1. Files already have the hammer and sensor calibration factor in them, but they are presented anyway just in case. % % % % % 2. G consists of 6 columns the first of which is the FRF that we want, so it will be extracted. % % % % % 3. f consists if the frequencies used in the G vector. % % % % % % % % [FILENAME,PATH] = uigetfile('*.txf'); % % % % FILENAME = [PATH FILENAME]; % % % % [f G] = Txfnew(FILENAME); % % % % % % % % % Accelerometer 1027 m/s^2/V for the small one. % % % % % 9 19.4 m/s^2/V for the slightly larger one % % % % % laser vibrometer 5e3 m/s/V % % % % % cap. probe = 2.5e5 m/V; % % % % % % % % % Hammer 925.93 N/V for the medium one with vinyl tip % % % % % 892.86 N/V for the medium one with steel tip % % % % % 46.9484 N/V for small hammer 086D80 serial# 19788 % % % % % % % % index = find(f>150); % % % % G = G(index,1); % % % % f = f(index); % % % % % % % % fmeas = f; % % % % FRFy = G; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Variables % %%%%%%%%%%%%% %%% ***Please set cutting coefficient order from low to high. There should be 3 values.*** % Specific cutting energy values Ktc_vec = 2395e6;%/700e6*[550e6 700e6 850e6]; % tangential cutting force coefficient, N/m^2 Krc_vec = 0.3*Ktc_vec; % radial cutting force coefficient, N/m^2 Kte_vec = 0*ones(1,length(Ktc_vec)); % tangential edge constant, N/m Kre_vec = 0*ones(1,length(Ktc_vec)); % radial edge constant, N/m

PAGE 114

114 % Tool description N = 4; % teeth, integer beta = 45; % helix angle, deg d = 10e3; % teeth diameter, m % Stability lobes Num_lobes = 14; ss = 5000:10:50000; % spindle speed, rpm % Cost info W = 100e3; % critical dimension of the cube part (formualtion may change) in m TLC = 1.9549e6; exp1 = 1.6265; % v exponent for TTL (exponents must be negative) exp2 = 0.1024; % ft exponent for TTL (exponents must be negative) exp3 = 0.2837; % b exponent for TTL (exponents must be negative) rm = 1; % machining cost per min ($/min) tch = 4/60; % tool changing time (min) Ct = 114; % Cost of a tool ($) R = 804.48; % Revenue per part % Machining specifications for test cases RDOC = [2 2.5 3 3.5 4 4.5 5]*1e3; b_vec = [1 1.5 2 2.5 3 3.5 4]*1e3; ft = [.05 .07 .09 .11 .15 .17 .19]*1e3; SS_vec = linspace(23000,39000,7); num_R = length(RDOC); num_b = length(b_vec); num_f = length(ft); num_S = length(SS_vec); num_tot = num_R*num_b*num_f*num_S; % RDOC = [2 3.5 5]*1e3; % Radial Depth of Cut in m % b_vec = [1 2.5 5]*1e3; % axial depth of cut, m % ft = [0.01 0.075 0.15]*1e3; % feed/tooth, m % SS_vec = [10000 24700 37000]; % Spindle Speeds in RPM rho = RDOC/d; % Radial Immersion as a fractional quantity not a percent % Milling direction used to determine entry and exit angle. %1 = upmilling and 2 = downmilling dir = 2; ADM = 0; % Axial depth where SLE is measured, m (note that 0 is the tip of the tool) % SLE Constraint

PAGE 115

115 SLE_MAX = 50; % micrometers Summary_mat = zeros(num_tot,7,length(Ktc_vec)); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Calling programs and other calculations for BIG_CNT = 1:length(Ktc_vec) % "MEAN" FMC data calculations Ktc = Ktc_vec(BIG_CNT); Krc = Krc_vec(BIG_CNT); Kte = Kte_vec(BIG_CNT); Kre = Kre_vec(BIG_CNT); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % DETERMINING ENTRY AND EXIT ANGLE % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if rho <= 0.5 % Entry and Exit angles (radians) if dir == 1 phi_in = zeros(1,length(rho)); % entry angle (upmilling) phi_out = acos(1rho*2); % exit angle (upmilling) elseif dir == 2 phi_in = piacos(1rho*2); % entry angle (downmilling) phi_out = pi*ones(1,length(rho)); % exit angle (downmilling) else fprintf('The direction you have chosen is incorrect please run again and select "up = 1" or "down = 2" for the dir variable. \ n'); return end elseif rho > 0.5 & rho<1 if dir == 1 phi_in = zeros(1,length(rho)); phi_out = phi_in + pi/2 + acos(1(rho0.5)*2); elseif dir == 2 phi_out = pi*ones(1,length(rho)); phi_in = phi_out pi/2 acos(1(rho0.5)*2); else fprintf('The direction you have chosen is incorrect please run again and select "up = 1" or "down = 2" for the dir variable. \ n'); return end elseif rho == 1 phi_in = zeros(1,length(rho)); phi_out = pi*ones(1,ength(rho)); elseif rho > 1 || rho <=0; fprintf('Please review the radial immersion of the cut, the current value is illogical. \ n'); return end

PAGE 116

116 phistart = phi_in*180/pi; % starting angle, deg phiexit = phi_out*180/pi; % exit end, deg %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% quant1 = length(RDOC)*length(ft)*length(b_vec)*length(SS_vec); % # all combinations quant2 = length(ft)*length(b_vec)*length(SS_vec); % # combinations except RDOC quant3 = length(b_vec)*length(SS_vec); % # ADOC and SS combinations quant4 = length(RDOC)*length(ft); % Combinations of RDOC and ft DTR = pi/180; % conversion Degrees to radians toc %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Stability %%%%%%%%%%%%%%%%%% Stability_data_mat = zeros(quant1,5); Stability_data_mat(:,1) = 1:quant1; for cnt = 1:length(RDOC) index = 1+quant2*(cnt1); Stability_data_mat(index:(index+quant21),2) = [RDOC(cnt)*ones(quant2,1)]; end SS_lobes_vec_data = zeros(length(RDOC),4501); Blim_lobes_vec_data = zeros(length(RDOC),4501); for cnt5 = 1:length(RDOC) [Combination_mat,SS_final,Blim_final] = stability_alt_trim(Ktc, Krc, N, phistart(cnt5)*DTR, phiexit(cnt5)*DTR, ss, Num_lobes, b_vec*1e3, SS_vec, fmeas, FRFx, FRFy); for cnt6 = 1:length(ft) index1 = 1 + quant3*(cnt61)+quant2*(cnt51); index2 = quant3 + quant3*(cnt61)+quant2*(cnt51); Stability_data_mat(index1:index2,3:5) = [Combination_mat]; end SS_lobes_vec_data(cnt5,:) = SS_final; Blim_lobes_vec_data(cnt5,:) = Blim_final; end toc

PAGE 117

117 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % SLE %%%%%%%%%%%%%%%%%% counter = 0; SLE_data_mat = zeros(quant1,7); for cnt1 = 1:length(RDOC) for cnt2 = 1:length(ft) for cnt3 = 1:length(b_vec) for cnt4 = 1:length(SS_vec) counter = counter+1; [e] = sle_navy_DCcorr_axialselect(fmeas, FRFy, Ktc, Krc, Kte, Kre, N, beta, d, phistart(cnt1), phiexit(cnt1), SS_vec(cnt4), b_vec(cnt3), ft(cnt2), ADM); sle = e*1e6; if abs(real(sle))<=SLE_MAX SLE_mult = 1; elseif abs(real(sle))>SLE_MAX SLE_mult = 0; end SLE_data_mat(counter,:) = [counter RDOC(cnt1) ft(cnt2) b_vec(cnt3) SS_vec(cnt4) real(sle) SLE_mult]; end end end end toc %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Ra %%%%%%%%%%%%%%%%%% Ra_max = 1; Ra_mat = zeros(quant1,2); counter = 0; for cnt1 = 1:length(RDOC) for cnt2 = 1:length(ft) for cnt3 = 1:length(b_vec) for cnt4 = 1:length(SS_vec) Ra = 10^3*(ft(cnt2)^2/(32*(5ft(cnt2)*4/pi))); if abs(Ra) > abs(Ra_max) Ra_mult =0; elseif abs(Ra) <= abs(Ra_max)

PAGE 118

118 Ra_mult =1; end counter = counter+1; Ra_mat(counter,:) = [counter Ra_mult]; end end end end toc %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Cost and Tool life encoded Profit %%%%%%%%%%%%%%%%%% [Profit_data_mat] = Profit_data_Toollife_encoded(W*1e3,d*1e3,N,RDOC*1e3,ft*1e3,b_vec*1e3,SS_vec,T LC,exp1,exp2,exp3,rm,tch,Ct,R); toc %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Multipliers %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Final_mult = Stability_data_mat(:,end).*SLE_data_mat(:,end).*Ra_mat(:,end); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Summary so far (last two columns will be empty aka zeros) % The summary matrix is compossed of 8 columns (# RDOC ft b_mean SS P_mean P_upper b_upper) Summary_mat(:,:,BIG_CNT) = [SLE_data_mat(:,1:5) Profit_data_mat(:,end) Final_mult]; toc end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % MAX Profit Calculation if length(Ktc_vec) == 3 Profit_encoded = 0.25*Summary_mat(:,6,1).*Summary_mat(:,7,1) + 0.5*Summary_mat(:,6,2).*Summary_mat(:,7,2) + 0.25*Summary_mat(:,6,3).*Summary_mat(:,7,3); Final_summary_mat = [Summary_mat(:,1:5,1) Profit_encoded]; elseif length(Ktc_vec)==1; Profit_encoded = Summary_mat(:,6,1).*Summary_mat(:,7,1) ; Final_summary_mat = [Summary_mat(:,1:5,1) Profit_encoded];

PAGE 119

119 end CHECK = Final_summary_mat; id_infeasible = find(Final_summary_mat(:,end) == 0); CHECK(id_infeasible,:) = []; id_max = find(CHECK(:,end) == max(CHECK(:,end))); MAX_Profit = max(CHECK(:,end)) MAX_Profit_parameters = CHECK(id_max,1:5) toc

PAGE 120

120 Tool Life % Ral E. Zapata % % Cost Contour Plots for machining optimization % % October 2006 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [Profit_data_mat] = Profit_data_Toollife_encoded(W,d,N,a,ft,b,RPM,TLC,exp1,exp2,exp3,rm,tch,Ct,R) % This function will provide the cost per part for all the desired data % points along with the input combinations for a, b, ft, and RPM that % caused them using a Taylor type tool life equation. % Note: Keeping track of units is essential!!!! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Calculations %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% v = RPM*pi*d/1000; % Surface velocity (m/min) counter = 0; Cost_data_mat = zeros(length(a)*length(b)*length(ft)*length(RPM),7); delta = 4.72; T2_enc = 30; T1_enc = T2_encdelta; T3_enc = T2_enc+delta; for cnt1 = 1:length(a) for cnt2 = 1:length(ft) for cnt3 = 1:length(b) for cnt4 = 1:length(v) counter = counter+1; L = W/b(cnt3)*(2*W/a(cnt1)*(W+d)+W); Lc = W/b(cnt3)*(W/a(cnt1)*(W+d)); f = RPM(cnt4)*N*ft(cnt2); T2 = TLC v(cnt4)^exp1 ft(cnt2)^exp2 b(cnt3)^exp3; % Tool life in minutes T2_ratio = T2/T2_enc; T1 = T1_enc*T2_ratio; T3 = T3_enc*T2_ratio; tm = L/f; % machining time in minutes tc = Lc/f; % Cutting time in minutes

PAGE 121

121 Cm1 = tm*rm + (tch*rm+Ct)*tc/T1; % Cost in Dollars Cm2 = tm*rm + (tch*rm+Ct)*tc/T2; Cm3 = tm*rm + (tch*rm+Ct)*tc/T3; P1 = RCm1; P2 = RCm2; P3 = RCm3; Pexp = P1*.25+P2*.5+P3*.25; Profit_data_mat(counter,:) = [counter a(cnt1) ft(cnt2) b(cnt3) RPM(cnt4) T1 T2 T3 P1 P2 P3 Pexp]; end end end end

PAGE 122

122 Fast Fourier Transform C ode Used In SLE Code % Computes the fft X of signal x and the corresponding frequency vector f given % the sampling frequency fs. % % [X,f]=spec(x,fs) % % [X,f]=spec(x,fs,'whole') returns values around the whole unit circle function [X,f]=spec(x,fs,whole) T=1/fs; N=length(x); %X=T*fft(x); X=fft(x); f=[0:fs/N:(11/(2*N))*fs]'; if nargin == 2 X=X(1:N/2+1,:); f=f(1:N/2+1,:); end

PAGE 123

123 APPENDIX C SUPER DIAGRAM CALCUL ATION CODES. The first of these codes presents the program used to compute and arrange the stability and surface location error information at each test point. This data is saved and then used in the second program that plot s the super diagram(s) according to the user s tolerance preferences.

PAGE 124

124 Data Gathering % Ral E. Zapata % Call program for super diagram data collection % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % This program calls the stability_alt_trim program to create stability % lobes and test a series of operating points to see if they are above or % below the stability limit (unstable or stable). % Afterwards the program inputs these same points into a surface location % error calculation to obtain the contours of levels of error. % These are then combined to form a diagram delimiting the areas that are % and are not of use according to the user selected bounds of surface % location error. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% close all clear all clc % Inputs. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Variables % %%%%%%%%%%%%% %%% ***Please set cutting coefficient order from low to high. There should be 3 values.*** % Specific cutting energy values % % Mean Real Values % Ktc_vec = 834e6; % tangential cutting force coefficient, N/m^2 % Krc_vec = 472e6; % radial cutting force coefficient, N/m^2 % Kte_vec = 2.62e3; % tangential edge constant, N/m % Kre_vec = 1.75e3; % radial edge constant, N/m % % Max Real Values 1mm cuts % Ktc_vec = 1040.6e6; % tangential cutting force coefficient, N/m^2 % Krc_vec = 618.7e6; % radial cutting force coefficient, N/m^2 % Kte_vec = 1.984e3; % tangential edge constant, N/m % Kre_vec = 1.604e3; % radial edge constant, N/m % Max Real Values 0.75mm cuts Ktc_vec = 700e6; % tangential cutting force coefficient, N/m^2 Krc_vec = 210e6; % radial cutting force coefficient, N/m^2 Kte_vec = 0; % tangential edge constant, N/m Kre_vec = 0; % radial edge constant, N/m

PAGE 125

125 % Tool description N = 4; % teeth, integer beta = 45; % helix angle, deg d = (1/4*25.4)*1e3; % teeth diameter, m % Stability lobes Num_lobes = 10; % number of stability lobes to be calculated. ss = 5000:10:50000; % spindle speed, rpm % Machining specifications for test cases ft = .075e3; % feed per tooth (m/tooth) RDOC = d/2; % Radial depth of cut. rho = RDOC/d; % Radial Immersion as a fractional quantity not a percent % Milling direction used to determine entry and exit angle. %1 = upmilling and 2 = downmilling dir = 1; ADM = 0; SLE_MAX = 30; %%%%%%%%%%%%%%%%%%%%%%% % FRF %%%%%%%%%%%%%%%%%%%%%%% % MODAL FIT (I created) n_teeth = 4; % # of teeth on the tool SS_best1 = 18000*2; % "best speed" arbitrary value for the peak around the first lobe. F_best1 = SS_best1/60*n_teeth; % % Tool description ky = [5e6]; % N/m zetay = [0.05]; % damping coefficient wny = [F_best1]*2*pi; % rad/s my = ky./(wny.^2); % kg cy = 2*zetay.*(my.*ky).^0.5; % Ns/m % Define ydirection FRF (could read in measurement from TXF) fmeas = 0:1:5000; % Hz w = fmeas*2*pi; % rad/s FRFy = (wny(1)^2/ky(1))./(wny(1)^2 w.^2 + i*2*zetay(1)*wny(1).*w); % this for loop adds up any other modes of vibration to the FRF if they % exist in the formulation above. for cnt = 2:length(ky) FRFy = FRFy + (wny(cnt)^2/ky(cnt))./(wny(cnt)^2 w.^2 + i*2*zetay(cnt)*wny(cnt).*w); end

PAGE 126

126 FRFx = FRFy; % for this test I am assumming the X and Y dynamics are the same. % figure % subplot(211) % plot(fmeas,real(FRFx),'b') % subplot(212) % plot(fmeas,imag(FRFx),'b') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % .TXF FILE CALL % Assumptions: % 1. Files already have the hammer and sensor calibration factor in them, but they are presented anyway just in case. % 2. G consists of 6 columns the first of which is the FRF that we want, so it will be extracted. % 3. f consists if the frequencies used in the G vector. % % Y direction % [FILENAME,PATH] = uigetfile('*.txf'); % FILENAME = [PATH FILENAME]; % [f G] = Txfnew(FILENAME); % % index = find(f>100 & f<6000); % G = G(index,1); % f = f(index); % % fmeas = f; % FRFy = G; % % % X direction % [FILENAME,PATH] = uigetfile('*.txf'); % FILENAME = [PATH FILENAME]; % [f G] = Txfnew(FILENAME); % % index = find(f>100 & f<6000); % G = G(index,1); % f = f(index); % % fmeas = f; % FRFx = G; % % figure % subplot(211) % plot(fmeas,real(FRFx)) % ylabel('Real (m/N)') % subplot(212) % plot(fmeas,imag(FRFx)) % xlabel('Frequency (Hz)') % ylabel('Imaginary (m/N)') %%%%%%%%%%%%%%%%%%%%%%%

PAGE 127

127 % Test parameters %%%%%%%%%%%%%%%%%%%%%%% % Define test grid. spinspeed = linspace(15000,45000,300); ADOC = linspace(.0001,6,300); % Create test vectors count = 0; for cnt2 = 1:length(ADOC) for cnt1 = 1:length(spinspeed) count = count+1; SS_test(count) = spinspeed(cnt1); ADOC_test(count)= ADOC(cnt2); end end clear count % SS_test = [5000 5000 5000 7500 7500 7500 9000 15000 15000 18000 18000 18000 18000]; % test point spindle speeds % ADOC_test = [.5 1 2 .5 1 2 1 .5 2 .5 1 2 3]/2; % test point axial depths %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Calculations and procedures/functions %%%%%%%%%%%%%%%%%%%%%% % "MEAN" FMC data calculations Ktc = Ktc_vec; Krc = Krc_vec; Kte = Kte_vec; Kre = Kre_vec; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % DETERMINING ENTRY AND EXIT ANGLE % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if rho <= 0.5 % Entry and Exit angles (radians) if dir == 1 phi_in = zeros(1,length(rho)); % entry angle (upmilling) phi_out = acos(1rho*2); % exit angle (upmilling) elseif dir == 2 phi_in = piacos(1rho*2); % entry angle (downmilling) phi_out = pi*ones(1,length(rho)); % exit angle (downmilling) else

PAGE 128

128 fprintf('The direction you have chosen is incorrect please run again and select "up = 1" or "down = 2" for the dir variable. \ n'); return end elseif rho > 0.5 & rho<1 if dir == 1 phi_in = zeros(1,length(rho)); phi_out = phi_in + pi/2 + acos(1(rho0.5)*2); elseif dir == 2 phi_out = pi*ones(1,length(rho)); phi_in = phi_out pi/2 acos(1(rho0.5)*2); else fprintf('The direction you have chosen is incorrect please run again and select "up = 1" or "down = 2" for the dir variable. \ n'); return end elseif rho == 1 phi_in = zeros(1,length(rho)); phi_out = pi*ones(1,ength(rho)); elseif rho > 1 || rho <=0; fprintf('Please review the radial immersion of the cut, the current value is illogical. \ n'); return end phistart = phi_in*180/pi; % starting angle, deg phiexit = phi_out*180/pi; % exit end, deg %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DTR = pi/180; % conversion Degrees to radians %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Stability %%%%%%%%%%%%%%%%% % This should output the stability conditions and operating conditions % selected for testing in a table or matrix form, along with the vectors % required to create a stability boundary plot (spindle speed and axial depth). tic [Combination_mat,SS_final,Blim_final] = stability_alt_trim(Ktc, Krc, N, phistart*DTR, phiexit*DTR, ss, Num_lobes,ADOC_test,SS_test, fmeas, FRFx, FRFy); figure plot(SS_final,Blim_final,'b');%,SS_test,ADOC_test,'r*'); xlabel('Spindle speed (rpm)') ylabel('Axial depth (mm)') title(num2str(RDOC)) axis([Inf Inf 0 10]) stab = Combination_mat(:,3);

PAGE 129

129 toc %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % SLE %%%%%%%%%%%%%%%%% quant1 = length(stab); SLE_data_mat = zeros(quant1,4); counter = 0; tic for cnt = 1:length(ADOC) [e] = sle_SSvec_capable(fmeas, FRFy, Ktc, Krc, Kte, Kre, N, beta, d, phistart, phiexit, spinspeed, ADOC(cnt)*1e3, ft, ADM); sle = e*1e6; for cnt1 = 1:length(sle) counter = counter+1; if abs(real(sle(cnt1)))<=SLE_MAX %acceptable SLE_mult = 0; elseif abs(real(sle(cnt1)))>SLE_MAX %unacceptable SLE_mult = 1; end SLE_data_mat(counter,:) = [spinspeed(cnt1) ADOC_test(cnt) real(sle(cnt1)) SLE_mult]; end end SLE = SLE_data_mat(:,4); toc %%%%%%%%%%%%%%%%%%%%%%%%% % Data for Diagram %%%%%%%%%%%%%%%%%%%%%%%%% Total = stab+SLE; Total_mat = zeros(length(ADOC),length(spinspeed)); for cnt = 1:length(spinspeed) Total_mat(:,cnt) = Total(1+(cnt1)*length(ADOC):cnt*length(ADOC)).'; end figure contourf(spinspeed,ADOC,Total_mat) xlabel('Spindle speed (rpm)') ylabel('Axial Depth (mm)') save stab_SLE_superdiagram_data_p25inchmill_fpt075_2.mat SLE_data_mat Combination_mat spinspeed ADOC

PAGE 130

130 Data Organization And Plotting close all clear all clc % Load previous stability and SLE calculation % load stab_SLE_superdiagram_data.mat % load stab_SLE_superdiagram_data_ALUMINUM.mat % load stab_SLE_superdiagram_data_ALUMINUM2.mat % load stab_SLE_superdiagram_data_ALUMINUM3.mat % load stab_SLE_superdiagram_data_ALUMINUM3_upmill1.mat % load stab_SLE_superdiagram_data_REALCOEFF_maxvalues1mmADOC.mat % load stab_SLE_superdiagram_data_REALCOEFF_maxvalues_p75mmADOC.mat % load stab_SLE_superdiagram_data_REALCOEFF_HIDEF_maxvalues_p75mmADOC.mat % For numerical examples in ASPE paper % Erroneous % load stab_SLE_superdiagram_data_p125inchmill_fpt15.mat % load stab_SLE_superdiagram_data_p125inchmill_fpt075.mat % Correct % load stab_SLE_superdiagram_data_p25inchmill_fpt15.mat load stab_SLE_superdiagram_data_p25inchmill_fpt075.mat % load stab_sle_data.mat % spinspeed = linspace(5000,20000,300); % ADOC = linspace(.0001,3,300); % recalculate SLE boundaries for different levels. (First was set to 50) % SLE_MAX = linspace(20,120,11); % SLE_MAX = [6.35 12.7 20 30 40 50]; % SLE_MAX = [30 40 50 60 70 80 90 ]; SLE_MAX = [30 70]; % SLE_MAX = [42]; SLE_data_mat_expanded = zeros(length(SLE_data_mat(:,1)),3+length(SLE_MAX)); SLE_data_mat_expanded(:,1:3) = SLE_data_mat(:,1:3); for cnt = 1:length(SLE_MAX) for cnt2 = 1:length(SLE_data_mat(:,1)) if abs(real(SLE_data_mat(cnt2,3)))<=SLE_MAX(cnt) %acceptable SLE_data_mat_expanded(cnt2,cnt+3) = 0; elseif abs(real(SLE_data_mat(cnt2,3) ))>SLE_MAX(cnt) %unacceptable SLE_data_mat_expanded(cnt2,cnt+3) = 1; end end

PAGE 131

131 end %%%%%%%%%%%%%%%%%%%%%%%%% % Data for Diagrams %%%%%%%%%%%%%%%%%%%%%%%%% % First for loop only makes sense if the SLE predictions were valid in the % unstable space. However, SLE predictions are based on the assumption that % the cut was stable, therefore all SLE data in the unstable area is % invalid. % for cnt51 = 1:length(SLE_MAX) % Total = 2*(1Combination_mat(:,3))+SLE_data_mat_expanded(:,3+cnt51); % % Total_mat = zeros(length(ADOC),length(spinspeed)); % % for cnt = 1:length(spinspeed) % Total_mat(:,cnt) = Total(1+(cnt1)*length(ADOC):cnt*length(ADOC)).'; % end % % figure % contourf(spinspeed,ADOC,3Total_mat) % xlabel('Spindle speed (rpm)') % ylabel('Axial Depth (mm)') % title(num2str(SLE_MAX(cnt51))) % colormap gray % end % Fixed conditions Total = zeros(1,length(Combination_mat(:,3))); for cnt51 = 1:length(SLE_MAX) for cnt52 = 1:length(Combination_mat(:,3)) % note that this may change depending on the input file that I chose. % The first trial this is ~= but after "fixing" the other program this % should be == for the first condition. if Combination_mat(cnt52,3) == 0 if SLE_data_mat_expanded(cnt52,3+cnt51) == 0 Total(cnt52) = 0; elseif SLE_data_mat_expanded(cnt52,3+cnt51) ~= 0 Total(cnt52) = 1; end elseif Combination_mat(cnt52,3) ~= 0 Total(cnt52) = 2; end end for cnt53 = 1:length(spinspeed) Total_mat(:,cnt53) = Total(1+(cnt531)*length(ADOC):cnt53*length(ADOC)).'; end

PAGE 132

132 figure contourf(spinspeed,ADOC,3Total_mat) xlabel('Spindle speed (rpm)') ylabel('Axial depth (mm)') title(num2str(SLE_MAX(cnt51))) colormap gray end

PAGE 133

133 APPENDIX D STABILITY UPDATING The Matlab code presented in this appendix completes the numerical example for stability updating provide d in Chapter 4. It accepts the inputs for stability as described and commented in the first section. Then, it generates the reference stability lobe an d the prior distribution. After generating the starting points it find a suitable optimal test point and proceeds to perform the prescribed number of test points required by the user.

PAGE 134

134 % function new_update2() close all clear all clc global switcher; global Blim_final; % global Blim_value; %%This code is set up to do the following. First, you input the initially %%known test points in test_grid below. This defines the initial prior %%distribution. Each time test_count clicks up one (each time through the first %%while loop), the code runs a mock updating for every possible spindle %%speed/axial depth combination for both stable and unstable results (200 x %%150 x 2 total mock updates). For each mock update, the resultant best %%place to actually perform cutting is calculated as well as the profit made by cutting there(to reduce complexity, it %%is assumed that we only will cut at a point which is known to be stable). % Then, the expected profit for each pair of mock tests (200 x 150 pairs) % is calculated (Expected Profit = Prior Prob of % Stability*(Resultant Profit|test_stable) +(1Prior Prob of % stability)*(Resultant Profit|test_unstable)). Then, the maximum expected profit is taken % over these testing possibilities (200 x 150) and the spindlespeed axial % depth which achieves this maximum is chosen for testing. Testing is then % performed (via simulation against the reference stability lobe). This % process is then repeated for the chosen number of tests (14 as the code % is currently written). tic switcher=0; b_min=0; b_max=1.2; b_number=200; b_space=b_min:(b_maxb_min)/b_number:b_max; l=b_maxb_min; s_min=5000; s_max=20000; s_number=150; s_space=s_min:(s_maxs_min)/s_number:s_max; k_const=.00005; sinterms=200; m_space=1:sinterms; number_of_tests=14; test_grid=sortrows([s_min b_min b_max; s_max b_min b_max],1); % Each row is a tested spindle speed, first column is that spindle speed, % second collumn is the highest stable axial depth tested (0 if no stable %tests at this spindle speed), third collumn is the lowest unstable spindle %speed tested (b_max if no stable tests at this spindle speed). %make sure to keep a row for the lowest considered spindle speed and the %highest considered spindle speed

PAGE 135

135 test_grid_indices(:,1)=(test_grid(:,1)s_min)*s_number/(s_maxs_min)+1; %%test grid converted into the indices test_grid_indices(:,2:3)=(test_grid(:,2:3)b_min)*b_number/(b_maxb_min)+1; spindle_speeds_tested=max(size(test_grid(:,1))); %%number of spindle speeds tested V_op=0; Value_of_K(1:3)=0; Ktc_vec = 2500e6; % tangential cutting force coefficient, N/m^2 Krc_vec = 0.3*Ktc_vec; % radial cutting force coefficient, N/m^2 Kte_vec = 0*ones(1,length(Ktc_vec)); % tangential edge constant, N/m Kre_vec = 0*ones(1,length(Ktc_vec)); % radial edge constant, N/m % Tool description N = 4; % teeth, integer beta = 45; % helix angle, deg d = 10e3; % teeth diameter, m % Stability lobes Num_lobes = 10; % number of stability lobes to be calculated. ss = s_space; % spindle speed, rpm % Machining specifications for test cases RDOC = d/2; % Radial depth of cut. rho = RDOC/d; % Radial Immersion as a fractional quantity not a percent % Milling direction used to determine entry and exit angle. %1 = upmilling and 2 = downmilling dir = 2; %%%%%%%%%%%%%%%%%%%%%%% % FRF %%%%%%%%%%%%%%%%%%%%%%% % MODAL FIT (I created) n_teeth = 4; % # of teeth on the tool SS_best1 = 18000*2; % "best speed" arbitrary value for the peak around the first lobe. F_best1 = SS_best1/60*n_teeth; % % Tool description ky = [5e6]; % N/m zetay = [0.05]; % damping coefficient wny = [F_best1]*2*pi; % rad/s my = ky./(wny.^2); % kg cy = 2*zetay.*(my.*ky).^0.5; % Ns/m % Define ydirection FRF (could read in measurement from TXF) fmeas = 0:0.1:5000; % Hz w = fmeas*2*pi; % rad/s

PAGE 136

136 FRFy = (wny(1)^2/ky(1))./(wny(1)^2 w.^2 + i*2*zetay(1)*wny(1).*w); % this for loop adds up any other modes of vibration to the FRF if they % exist in the formulation above. for cnt = 2:length(ky) FRFy = FRFy + (wny(cnt)^2/ky(cnt))./(wny(cnt)^2 w.^2 + i*2*zetay(cnt)*wny(cnt).*w); end FRFx = FRFy; % for this test I am assumming the X and Y dynamics are the same. % % figure % % subplot(211) % % plot(fmeas,real(FRFx),'b') % % subplot(212) % % plot(fmeas,imag(FRFx),'b') %%%%%%%%%%%%%%%%%%%%%%% % Test parameters %%%%%%%%%%%%%%%%%%%%%%% % Define test grid. spinspeed = 10000; ADOC = [.2 ]; % Create test vectors count = 0; for cnt1 = 1:length(spinspeed) for cnt2 = 1:length(ADOC) count = count+1; SS_test(count) = spinspeed(cnt1); ADOC_test(count)= ADOC(cnt2); end end % SS_test = [5000 5000 5000 7500 7500 7500 9000 15000 15000 18000 18000 18000 18000]; % test point spindle speeds % ADOC_test = [.5 1 2 .5 1 2 1 .5 2 .5 1 2 3]/2; % test point axial depths %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Calculations and procedures/functions %%%%%%%%%%%%%%%%%%%%%% % "MEAN" FMC data calculations Ktc = Ktc_vec; Krc = Krc_vec; Kte = Kte_vec; Kre = Kre_vec; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

PAGE 137

137 % DETERMINING ENTRY AND EXIT ANGLE % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if rho <= 0.5 % Entry and Exit angles (radians) if dir == 1 phi_in = zeros(1,length(rho)); % entry angle (upmilling) phi_out = acos(1rho*2); % exit angle (upmilling) elseif dir == 2 phi_in = piacos(1rho*2); % entry angle (downmilling) phi_out = pi*ones(1,length(rho)); % exit angle (downmilling) else fprintf('The direction you have chosen is incorrect please run again and select "up = 1" or "down = 2" for the dir variable. \ n'); return end elseif rho > 0.5 & rho<1 if dir == 1 phi_in = zeros(1,length(rho)); phi_out = phi_in + pi/2 + acos(1(rho0.5)*2); elseif dir == 2 phi_out = pi*ones(1,length(rho)); phi_in = phi_out pi/2 acos(1(rho0.5)*2); else fprintf('The direction you have chosen is incorrect please run again and select "up = 1" or "down = 2" for the dir variable. \ n'); return end elseif rho == 1 phi_in = zeros(1,length(rho)); phi_out = pi*ones(1,ength(rho)); elseif rho > 1 || rho <=0; fprintf('Please review the radial immersion of the cut, the current value is illogical. \ n'); return end phistart = phi_in*180/pi; % starting angle, deg phiexit = phi_out*180/pi; % exit end, deg %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DTR = pi/180; % conversion Degrees to radians clear i [Combination_mat,SS_final,Blim_final] = stability_alt_trim(Ktc, Krc, N, phistart*DTR, phiexit*DTR, ss, Num_lobes,ADOC_test,SS_test, fmeas, FRFx, FRFy); toc % Initialize variables to save time

PAGE 138

138 phi = zeros(test_grid_indices(2,1),b_number+1); phi_hat = phi; b = zeros(1,sinterms); b_hat = b; CDF_sizing = zeros(s_number+1,length(phi(1,:))); init_test_grid_size=size(test_grid,1); endpoint = max(init_test_grid_size+1:init_test_grid_size+number_of_tests); CDF_hist = zeros(length(phi(1,:)),s_number+1,endpoint2); %%this is the code for calculating the updated probability distributions for test_count=init_test_grid_size+1:init_test_grid_size+number_of_tests clear test_grid_indices test_grid_indices(:,1)=(test_grid(:,1)s_min)*s_number/(s_maxs_min)+1; test_grid_indices(:,2:3)=fix((test_grid(:,2:3)b_min)*b_number/(b_maxb_min)+1); spindle_speeds_tested=max(size(test_grid(:,1))); for j=1:test_grid_indices(2,1) for k=1:b_number+1; phi(j,k)=sin(pi*b_space(k)/l); end phi(j,:)=phi(j,:)/trapz(b_space,phi(j,:)); end %%build the rightward traveling Brownian motion... i2 counts the number of %%tested spindle speeds crossed (the space between two tested spindle %%speeds is calculated as a group) for i2=3:size(test_grid,1) j=test_grid_indices(i21,1); for m=1:sinterms; b(m)=2/l*trapz(b_space(test_grid_indices(i21,2):test_grid_indices(i21,3)),phi(j,test_grid_indices(i21,2):test_grid_indices(i21,3)).*sin(m*pi*(b_space(test_grid_indices(i21,2):test_grid_indices(i21,3))b_min)/l)); %%Calculating Fourier series coefficients to the diffusion equation/Brownian %%Motion. sinterms is the number of terms in the series expansion, %%b(m) is the fourier coefficients. end j=j1; while(j
PAGE 139

139 %%summing the series expansion end phi(j,:)=phi(j,:)/trapz(b_space,phi(j,:)); end end %%build phi_hat_1(t,x) j=test_grid_indices(spindle_speeds_tested,1)+1; while(j>test_grid_indices(spindle_speeds_tested1,1)) j=j1; for k=1:b_number+1 phi_hat(j,k)=sin(pi*b_space(k)/l); end phi_hat(j,:)=phi_hat(j,:)/trapz(b_space,phi_hat(j,:)); end %%build the leftward traveling Brownian Motion. This mirrors the above %%code except now we start at the highest spindle speed, i2 counts the %%number of tested spindle speeds we've crossed starting from the right and %%traveling left. See the explanation for the rightward traveling Brownian %%motion for i2=1:size(test_grid,1)2 j=test_grid_indices(spindle_speeds_testedi2,1); for m=1:sinterms b_hat(m)=2/l*trapz(b_space(test_grid_indices(spindle_speeds_testedi2,2):test_grid_indices(spindle_speeds_testedi2,3)),phi(j,test_grid_indices(spindle_speeds_testedi2,2):test_grid_indices(spindle_speeds_testedi2,3)).*sin(m*pi*(b_space(test_grid_indices(spindle_speeds_testedi2,2):test_grid_indices(spindle_speeds_testedi2,3))b_min)/l)); end j=j+1; while(j>test_grid_indices(spindle_speeds_testedi21,1)) j=j1; for k=1:b_number+1 for m=1:sinterms phi_hat_terms(j,k,m)=b_hat(m)*exp(m^2*pi^2*k_const*(s_space(j)+s_space(test_grid_indices(spindle_speeds_testedi2,1)))/l^2)*sin(m*pi*(b_space(k)b_min)/l); end phi_hat(j,k)=sum(phi_hat_terms(j,k,:)); end phi_hat(j,:)=phi_hat(j,:)/trapz(b_space,phi_hat(j,:)); end end %%build and normalize the probability density function, which equals %%(leftward traveling Brownian Motion)*(rightward traveling Brownian %%Motion). clear mu CDF % mu = CDF_sizing; % CDF = CDF_sizing; % % for j=1:s_number+1 % mu(j,:)=phi(j,:).*phi_hat(j,:); % mu(j,:)=mu(j,:)/trapz(b_space,mu(j,:)); %%Normalize the distribution

PAGE 140

140 % CDF(j,:)=cumtrapz(b_space,mu(j,:)); %%Integrate to calculate the cumulative distribution function, aka the probability of stability. % end mu = phi.*phi_hat; denominator = trapz(b_space,mu,2); for cnt = 1:length(denominator) mu(cnt,:) = mu(cnt,:)/denominator(cnt); end CDF=cumtrapz(b_space,mu,2); figure subplot(311) contourf(s_space,b_space,phi.') ylabel('Axial depth (mm)') subplot(312) contourf(s_space,b_space,phi_hat.') ylabel('Axial depth (mm)') subplot(313) contourf(s_space,b_space,mu.') ylabel('Axial depth (mm)') xlabel('Spindle speed (rpm)') %%toc %%This is the code for the profit calculations r_m=1; %%$/min W=500; %%mm C_t=114; %%$ t_ch=4/60; %%min N=4; %%number of teeth d=10; %%diameter of tool in mm a=5; %%radial depth of cut (mm) L=@(b) W/b*(2*W/a*(W+d)+W); %%a is radial depth (mm), b is axial depth (mm), d is tool diameter (mm), L is length of the workpiece L_c=@(b) W/b*(W/a*(W+d)); %%mm, L_c is cutting length t_m=@(b,f) L(b)/f; %%f is feed rate (mm/min), t_m is in machining time (while cutting and while just moving tool) in minutes t_c=@(b,f) L_c(b)/f; %%t_c is cutting time in minutes f_t=@(f,s) f/s/N; %%s is in rpm f_t is mm/tooth (feed per tooth) v=@(s)pi*d*s/1000; %% cutting speed in m/min T=@(b,f,s) 19.9549*10^6*v(s)^1.6265*f_t(f,s)^.1024*b^.2837; %% tool life, v is m/min, ft is mm/tooth b in mm C_m=@(b,f,s) t_m(b,f)*r_m+(t_ch*r_m+C_t)*t_c(b,f)/T(b,f,s); %%cost of machining ADOC=b_space; %% axial depth of cut (mm) SS=s_space; %% spindle speed (rpm) f=10000; Cutting_speed=v(SS);

PAGE 141

141 C_m_space = zeros(size(SS,2),size(ADOC,2)); for i2=1:size(SS,2) for j=1:size(ADOC,2) C_m_space(i2,j)=C_m(ADOC(j),f,SS(i2)); %%calculates the cost of cutting at each spindle speed/axial depth combination after each test is completed, (you will have to look at the lines 298312 above to figure out how this was calculated) end end Profit=max(max(C_m_space))C_m_space; %%the maximum achievable profit Utility=Profit.*(1CDF)+CDF*max(Value_of_K); %%Value_of_K is a vector of how much profit is achieved at each previously tested point (since only cutting conditions known to be stable are considered, max(Value_of_K) is the profit if the new test is unstable) V_op=max(max(Utility)); ADOC_op_index=floor(find(Utility==V_op)/size(s_space,2))+1; SS_op_index=find(Utility==V_op)(ADOC_op_index1)*size(s_space,2); if(SS_op_index==0) SS_op_index=size(s_space,2); ADOC_op_index=ADOC_op_index1; end ADOC_op=b_space(ADOC_op_index); SS_op=s_space(SS_op_index); if(max(SS_op==test_grid(2:size(test_grid,1)1,1))==0) test_grid(size(test_grid,1)+1,1)=SS_op; test_grid(size(test_grid,1),2)=b_min; test_grid(size(test_grid,1),3)=b_max; elseif(test_count==3) test_grid(size(test_grid,1)+1,1)=SS_op; test_grid(size(test_grid,1),2)=b_min; test_grid(size(test_grid,1),3)=b_max; else test_grid(find(SS_op==test_grid(2:size(test_grid,1)1,1))+1,1)=SS_op; end %%updating the test_grid ADOC_test=ADOC_op; SS_test=SS_op; switcher=1; Blim_final(SS_op_index); if(ADOC_test<=Blim_final(SS_op_index)) stab=1; else stab=0; end if(max(SS_op==test_grid(2:size(test_grid,1)1,1))==0) if(stab==1) test_grid(size(test_grid,1),2)=ADOC_op; else test_grid(size(test_grid,1),3)=ADOC_op; end else if(stab==1) test_grid(find(SS_test==test_grid(2:size(test_grid,1)1,1))+1,2)=ADOC_op;

PAGE 142

142 else test_grid(find(SS_test==test_grid(2:size(test_grid,1)1,1))+1,3)=ADOC_op; end end test_grid=sortrows(test_grid,1); Value_of_K(test_count)=V_op; CDF_hist(:,:,test_count2) = CDF'; end toc figure surf(SS,ADOC,CDF','LineStyle', 'none') hold on plot3(SS,Blim_final,ones(length(SS),1),'Color', 'w') % CONT = linspace(0,1,100); % figure % view(2) % subplot(221) % view(2) % surf(SS,ADOC,1CDF_hist(:,:,1),'LineStyle','none') % hold on % plot3(SS,Blim_final,ones(length(SS),1),'Color','w') % axis([s_min s_max b_min b_max 0 1]) % view(2) % subplot(222) % view(2) % surf(SS,ADOC,1CDF_hist(:,:,5) ,'LineStyle','none') % hold on % plot3(SS,Blim_final,ones(length(SS),1),'Color','w') % axis([s_min s_max b_min b_max 0 1]) % view(2) % subplot(223) % view(2) % surf(SS,ADOC,1CDF_hist(:,:,9),'LineStyle','none' ) % hold on % plot3(SS,Blim_final,ones(length(SS),1),'Color','w') % view(2) % axis([s_min s_max b_min b_max 0 1]) % view(2) % subplot(224) % view(2) % surf(SS,ADOC,1CDF_hist(:,:,13),'LineStyle','none') % hold on % plot3(SS,Blim_final,ones(length(SS),1),'Color','w') % axis([s_min s_max b_min b_max 0 1]) % view(2) toc

PAGE 143

143 figure subplot(221) contourf(SS,ADOC,1CDF_hist(:,:,1),CONT,'LineStyle', 'none') hold on; plot(SS,Blim_final,'Color', 'w') axis([s_min s_max b_min b_max]) title('A') ylabel('Axial depth (mm)') subplot(222) contourf(SS,ADOC,1CDF_hist(:,:,5),CONT,'LineStyle', 'none') hold on plot(SS,Blim_final,'Color', 'w') axis([s_min s_max b_min b_max]) title('B') subplot(223) contourf(SS,ADOC,1CDF_hist(:,:,9),CONT,'LineStyle', 'none' ) hold on plot(SS,Blim_final,'Color', 'w') axis([s_min s_max b_min b_max]) title('C') xlabel('Spindle speed (rpm)') ylabel('Axial depth (mm)') subplot(224) contourf(SS,ADOC,1CDF_hist(:,:,13),CONT,'LineStyle', 'none') hold on plot(SS,Blim_final,'Color', 'w') axis([s_min s_max b_min b_max]) title('D') xlabel('Spindle speed (rpm)')

PAGE 144

144 LIST OF REFERENCES [1] Sherrif S Hisyam M, Kurhiawan D, Ovady E Performance evaluation of vegetable o il as an alternative cutting l ubricant Trans NAMRI/ SME 2009; 37: 914. [2] Cui Y, Fussel B Jeranol R, Esterling J Tool wear monitoring for milling by t racking cutting force model coefficients Trans NAMRI/SME 2009; 37: 61320. [3] Marusich T, Usui S, Lankalapalli S Marusich K, Garaud S, Zamorano L. P rediction of distortion of t hinwalled machined components Trans NAMRI/SME 2009; 37: 22936. [4] Schmitz T, Ziegert J Canning J, and Zapata R. Case S tudy: A c omparison of e rror s ources in milling Prec Eng 2008;32( 2): 126 133. [5] Schmitz T, Duncan G. T hree c omponent r eceptance c oupling s ubstructure a nalysis for t ool point dynamics prediction JMSE 2005; 127( 4): 781 90. [6] Altintas Y. Manufacturing a utomation. Cambridge, UK: Cambridge University Press; 2000. [7] Arnold R. The mechanism of t ool vibration in the cutting of steel In: Proceedings of the Institution of Mechanical Engineers 1946; 154( 4) : 26184. [8] Tobias S, Fishwick W. T he c hatter of l athe t ools under orthogonal c utting c onditions Trans ASME 1958 ; 80: 107988 [9] Tobias S, Fishwick W. T heory of r egenerative m achine t ool chatter The Engineer London 1958; 205:16 23. [10] Tobias S. M achine tool vibration Glasgow, Scotland: Blackie and Sons Ltd; 1965. [11] Mer rit H. Theory of self excited machine t ool chatter J of Eng for Ind 1965; 87( 4): 447 54. [12] Insperger T, Mann B, Stpn G and Bayly, P S tability of upmilling and downmilling Part 1: Alternative Analytical Methods IJMTM 2003; 43: 25 34. [13] Tlusty J Polocek M. The stability of the m achine t ool a gainst s elf excited vibration in m achining. In: Proceedings of the International Research in Production Engineering Conference Pittsburgh, PA 1963: 46574. [14] Altintas Y Budak E. A nalytical prediction of s tability l obes in m illing Annals of the CIRP 1995; 44( 1) : 35762. [15] Budak E and Altintas Y. A nalytical prediction of chatter stability conditions for multidegree of freedom systems in milling, P art I : M odeling, P art II: A pplications JDSMC, Trans of the ASME 1998 ; 120: 2236. [16] Ewins, D, Gleeson, P. A m ethod for m odal i dentification of l owly damped structures Journa l of Sound and Vibration 1982; 84( 1): 5779

PAGE 145

145 [17] Mann B, Young K, Schmitz T Bartow M, Bayly P. Machining accuracy due to tool or workpiece vibrations In: Proceedings of IMECE Washington, D.C., 2003 : 18. [18] Mann B, Bayly P, Davies M, Halley J. Limit c ycles, bifurcations and a ccuracy of the m illing process. J of Sound and Vibration 2004 ; 277 : 3148. [19] Schmitz, T. and Mann, B. Closed f orm s olutions for s urface l ocation e rror in m illing IJMTM, 2006; 46: 136977 [20] Alauddin M, El Baradie M and Hashmi M. Computer a ided a nalysis of a s urface r oughness m odel for end m illing JMPT 1995; 55: 123127. [21] Tlusty J. Manufacturing processes and equipment. Upper Saddle River, NJ: Prentice Hall; 2000. [22] Taylor F On the ar t of c utting metal Trans of the ASME 1907; 28: 31248 [23] Stephenson D, Agapiou J. Metal cutting theory and practice New York City, NY: Ma rcel Dekker; 1997. [24] H oward R, Matheson J. Influence diagrams. Decision Analysis 2005; 2( 3 ): 12743. [25] Spetzl er C, von Holstein C Probability e ncoding in decision a nalysis Management Science 1975; 22: 34058 [26] Smith J. Moment m ethods for decision a nalysis Management Science 1993; 39( 3): 340 58. [27] McNamee P Celona J. Decision analysis for the professional 3rd ed. Menlo Park, CA: SmartOrg Inc.; 2001. [28] Abbas A Entropy methods for univariate distributions in decision analysis In: 23rd International Workshop on Bayesian Inference and Maximum Entropy Methods in Science and Engineering, Moscow; 2002 (on CD) [29] Abbas A Maximum entropy distributions between upper and lower bounds In: 25th International Workshop on Bayesian Inference and Maximum Entropy Methods in Science and Engineering, San Jose, CA; 2005 (on CD) [30] Abbas A. Entropy methods for joint distributions in decision analysis. IEEE Trans on Engineering Management 2006; 53( 1):14659. [31] Ermer D, Patel D Maximization of the production rate with constraints by linear programming and sensitivity analysis In: Proceedings of NAMRC 1974; 2: 43649. [32] Ermer, D.S. Optimization of constrained machining economics pr oblem by geometric programming. ASME J of Engineering for Industry 1971; 93 : 106772. [33] Hati S, Rao S. Determination of optimum machining conditions determinist ic and probabilistic approaches ASME J of Engineering for Industry 1976; 98(1) : 3549.

PAGE 146

146 [34] Hitomi K. Analysis of optimal machining spee ds for automatic manufacturing. IJPR 1989; 27: 168591. [35] Shalaby M, Riad M. A linear optimization model for s ingle pass turning operations In: Proceedings of the 27th International MATADOR Conference 1988;27:2 3135. [36] Armarego E, Smith A Wang J. Constrained optimization strategies and cam soft ware for single pass peripheral milling IJPR 1993; 31( 9):213960. [37] Challa K, Berra P. Automated planning and optimization of machining processes: a systems approach Computers and Industrial Engineering 1976; 1: 3545. [38] Ju an H, Yu S, Lee B. The optimal cutting parameter selection of production cos t in HSM for SKD61 tool steels IJMTM 2003; 43( 7):679 86 [39] Iwata K Murostu, Y, I watsubo T Fujii S. A probabilistic approach to the determination of the optimum cutting. ASME J of Engineering for Industry 1972; 94: 1099107 [40] Philips D Beightler C. Optimization in tool engineering using geometric programming. IIE Trans 1970;2( 4) : 355 60. [41] Gopala krishnan B, Al Khayyal F. Machine parameter selection for turning with constraints: an analytical approach based on geometric programming. IJPR 1991; 29(9) : 1897908. [42] W alvekar A, Lambert B An Application of geometric programming to machining variable selection IJPR 1970; 8(3) : 24145. [43] Tandon V, E l Mounayri H, Kishawy H. NC end milling optimization using evolutionary computation. IJMTM 2002; 42( 5): 595 605. [44] Sonmez A, Baykasoglu A, Turkay D, Filiz I. Dynamic optimization of multi pass milling operations via geometric programming. IJMTM 1999; 39( 2):297320. [45] Wang J, Armarego E. Computer aided optimization of multiple constr aint single pass face milling operations Mach Science and Tech 2001; 5( 1): 7799. [46] Armarego E Smith A, Wang J. Computer aided constrained optimization analyses and strategies for multipass helical tooth milling operations. Annals of the CIRP 1994; 43( 1):43742. [47] Boothroyd G Rusek P. Maximum rate of profit criteria in machining ASME J of Engineering for Industry 1976; 98: 217 20. [48] Taylor F. Tool life studies for machining economics ASME 1893. [49] Gilbert W. The economics of machining. In: Machining Theory and Practice ASM 1950:46585

PAGE 147

147 [50] Gil bert, W. Economics of machining, in Machining with Carbides and Oxides. NewYork City New York: ASTME/McGraw Hill;1962. [51] Brewer R. On the economics of the basic turning operation. ASME J of Engineering for Industry 1958; 80: 147989 [52] Brown R. On the selecti on of economic machining rates. IJPR 1962; 11: 1 22. [53] McCullough E. Economics of multitool lathe operations ASME J of Engineering for Industry, 1963; 85: 40214. [54] Ham I. Economics of Machining : analyzing optimum machining condition s by computers. SP64 60, ASTME: 1964. [55] Shaw M. Optimum selection o f machine tool speeds and feeds. IJMTD 1965; 5: 2534. [56] Wu S Ermer D Hill W. A n exploratory study of taylor's tool life equation by power transformation ASME J of Engineering for Industry, 1964; 86 : 8192. [57] Wu S, Ermer D. Maximum profit rate as the criterion in the determination of the optimum cutting conditions. ASME J of Engineering for Industry 1966; 88: 43542. [58] Ermer D. Engineering statistics ap plied in production engineering. In: N. Am. Meta lworking Res. Conf. I, Ontario, Canada : 1973; 3. [59] Meng Lim E, Menq C H Integrated p lanning for precision m achining of c omplex s urfaces, P art 1: cutting path and f eedrate optmization. IJMTM 1967; 37( 1):61 75. [60] Makhanov S, Batanov D, Bohez E, Sonthipaumpoon K, Anotaipaiboon W, Tabuca non M. On t ool path optimization of a m illing robot. Computers and Industrial Engineering 2002; 43: 455 72. [61] Toh, C.K. Design, e valuation and optimization of c utter path s trategies when high s peed m achining hardened m old and die materials Materials and Design 2006; 26: 51733. [62] Komanduri R, Merchant M Shaw M S ymposium on the U.S. contributions to machining and grinding research in the 20th century. ASME Book No. AMR 126 or Applied Mechanics Review" 1993; 46( 3) :70 132 [63] El Mounayri H, Dugla Z, Deng H. Prediction of surface roughness in milling using swarm intelligence IEEE 2003; 2207. [64] Akturk M Avci S. Tool a llocation and m achining c onditions optimization for CNC machines. European Journal of Operational Research 1996; 94 (2 ) : 335 48. [65] Peres C Guerra R, Haber R, Alique A, Ros S. Fuzzy m odel and hierarchical f uzzy c ontrol i ntegration: an a pproach for m illing p rocess optimization Computers in Industry 1999; 39(3) : 199207.

PAGE 148

148 [66] Wang Z, Rahman M, Wong Y, Sun J. Optimization of m ultipass m illing using parallel genetic a lgorithm and parallel genetic s imulated a nnealing. IJMTM 2005; 45(15) : 1726 34. [67] Wan g J. Computer a ided e conomic optimization of e ndm illing operations IJPE 1998; 54(3) : 30720. [68] Lin T. Optimization technique for face milling stainless steel with multip le performance characteristics. IJAMT 2002; 19(5) : 330 5. [69] Tolouei Rad M, Bidhendi I. On the optimization of machining par ameters for milling operations. IJMTM 1997; 37( 1): 116. [70] Tzeng Y. A hybrid a pproach to optimize m ultiple performance c haracteristics of high s peed c omputerized numerical c ontrol m illing t ool steels. Materials and Design 2007; 28(1 ) : 3646. [71] Kurdi M. Robust multicriteria optimization of surface location error and material removal rate in high speed milling under uncertainty. Ph.D. D issertation University of Florida : 2005. [72] Kim H Schmitz T. Bivariate uncertainty a nalysis for i mpact t esting Measurement Science and Technology 2007; 18: 356571. [73] Tsa i M, Lee B Yu S. A predicted m odeling of t ool l ife of high s peed milling for SKD 61 tool steel. IJAMT 2005; 26(7 8) : 71117. [74] Sandvik Coromant, Rotating tools and inserts handbook, LITCAT 04 R. [75] DeGarm o E, Black J Kohser R. Materials and processes in manufacturing 9th Ed. Hoboken, NJ : Wiley and Sons ; 2003. [76] Taylor B Kuyatt C. Guidelines for evaluating and expressing the uncertai nty of NIST measurement results. NIST Technical Note 1297, 1994 Edition. [77] Ross S. Introduction to probability models 9th Edition. Burlington, MA: Academic Press, Elsevier; 2007. [78] Nagasawa, M. Schrodinge r equations and diffusion theory. Basel, Switzerland: Birkhauser Verlag; 1993. [79] Kac M. On distributions of certain W iener functionals Trans of the American Mathematical Society 1949; 65( 1): 113. [80] Trent E, Wright P. Metal c utting 4th Ed., Boston, MA : Butterworth Heinemann; 2000. [81] Hanna N Tobias S. A theory of nonlinear regenerative chatter. ASME J of Engineering for Industry 1974; 96: 24755.

PAGE 149

149 [82] Insperger T, Stpn G. Semi discretization method for delayed systems. IJNME 2002; 55( 5):50318. [83] Cheng, C.H. I mproved prediction of spindle holder tool frequency response functions Ph.D. Diss ertation, University of Florida; 2007.

PAGE 150

150 BIOGRAPHICAL SKETCH The author was born in San German, Puerto Rico in 1981. After finishing his high school education in 1999, he pursued his undergraduate degree at the University of Puerto Rico Mayaguez Campus in his hometown of Mayaguez. After four and a half years of scho ol, research, internship experiences, and fun in the sun he obtained his bachelors degree in mechanical engineering with high honors. His internship experience provided him with various job opportunities, which he turned down to pursue his graduate educat ion at the University of Florida. After a tumultuous first year, the author joined the Machine Tool Research Center (MTRC) with his current advisor Dr. Tony L. Schmitz and completed his m asters degree in 2006. H e immediately began work on his Ph.D. After completing his Ph.D., the author will continue working and teaching in the fields of manufacturing and optimization.