The Effect of Variable Time Delay on Stability in Orthogonal Cutting

1

The Pennsylvania State University The Graduate School Department of Mechanical Engineering

THE EFFECT OF VARIABLE TIME DELAY ON STABILITY IN ORTHOGONAL CUTTING

A Thesis in Mechanical Engineering  by Thomas Bernard O'Reilly

Submitted in Partial Fulfillment of the Requirements for the Degree of 

Master of Science

MAY 1998

2

3

Chapter 1

OVERVIEW

1.1. Background

Many industries depend upon machining processes to make a wide variety of  components. For economic reasons manufacturers want to produce parts as rapidly rapidly as  possible while being sure that the parts meet the design specifications for dimensional tolerances and surface finish. finish. Machining processes involve a machine tool such as a milling machine or a lathe, which holds a cutting tool and a workpiece and moves these two parts relative to one another. Machining processes all involve removing material from the workpiece so the goal of producing parts as quickly as possible is met by maximizing the material removal rate, which can be done by increasing the size of the chip being removed or by increasing the speed of the machine. In addition to this goal, companies want to be able to operate their machines with minimal downtime for repairs or tool changes. Several factors stand in the way of these objectives. One problem is static deformation of the machine tool structure structure due to cutting forces. Static deformation can lead to form errors and result in systematic errors errors in the sizes of parts produced. This  problem can be avoided by designing machine tools to be very stiff. Machine tool

4 manufacturers typically try to produce tools too ls with a static stiffness close to 1,000,000 lbs/in. Actual loop stiffnesses (the (the stiffness of the machine tool measured from the cutting tool to the part of the workpiece being machined) tend to be lower, but this type of error is generally well understood and can be dealt with in a number of ways. More difficult to deal with are dynamic displacements of either the machine tool or the workpiece. In addition to static deflections caused by mean mean cutting forces there is always the possibility that there will will be some vibration in the system. system. Such vibrations can begin many different ways, for example the spindle may not be balanced perfectly or  the slideways may not be perfectly straight. straight. Imperfections in the workpiece, such as voids or inclusions can also start vibrations. The cause of the vibrations is not critical. What is important important is that oscillations  between the tool and workpiece will cause the tool to leave a wa vy surface behind on the workpiece. In most machining processes the tool will cut through at least least a portion of the surface it just created. If that surface has any waviness on it the thickness of the chip  being removed by the cutting tool will vary. This variation in chip load will cause variations in the cutting and thrust forces exerted on the tool by the workpiece, causing the tool to oscillate again. In some cases the magnitude of the oscillations will decay with each pass of the tool due to damping in the process or structure. In other cases, the amplitude of the the oscillations can grow with successive passes passes of the cutting tool. This is a type of  instability known as regenerative chatter, and it can lead to excessive tool wear and very  poor surface quality. Oscillations due to this phenomenon can even become large enough

5 to cause the tool to leave the workpiece during a cut. Severe chatter marks are difficult to remove from a workpiece once they have been created due to the regenerative effect. In some cases it can be more economical to scrap the part and begin again. As a result, prediction of what conditions cause regenerative chatter is an important topic in the machine tool industry. industry. With more and more more automated machining centers being put into use the detection of chatter has also become an important topic as there is a need for machine tools to automatically detect and avoid chatter by altering  process parameters such as spindle speed or depth of cut. The stability limit is affected by many factors, such as depth of cut, tool geometry, feed rate and the stiffness and damping of the machine tool itself. itself. Also a factor in chatter are the properties of the material being machined. Hard materials like steel, titanium or ceramics are more prone to chatter than softer materials such as aluminum. As a result, chatter is a major impediment to high-speed machining of  difficult materials.

1.2. Problem Statement

Models of regenerative chatter make many simplifications to the mechanics of the cutting process and the dynamics of the system system response. Although these models provide useful clues about the behavior of a machine tool system they may neglect important details about the characteristics of the system.

6 For example, the vast majority of cutting models, like the one in Figure 1.1, only consider the effect of tool motion motion perpendicular to the surface of the workpiece. This is a good approach to build a basic model as it is this motion that will have the greatest impact on the profile of the surface left behind and on the force exerted on the cutting tool. From an analytical standpoint this is also a very attractive approximation as it results in a system with a known, constant delay period. Such a system is much easier to analyze than one with a variable delay period.

Cutting force vector  Chip flow

 N , Spindle speed in RPM

Θ

 x k  x m x c x

Figure 1.1: Single degree-of-freedom turning model.

However, in a real machine relative motion of the tool will likely have some component parallel to the surface of the workpiece. In addition the spindle spindle speed can be deliberately varied during a cutting process to improve stability. stability. In either of these cases

7 the delay period will not be a constant length of time. It will in fact be dependent on the magnitude of the tool displacement parallel to the surface, the displacement of the tool  perpendicular to the surface the last time it passed over that spot on the workpiece and any variations in spindle spindle speed. Some research has been performed in the the area of  variable speed machining, but most work in this area is still done on a trial and error   basis. Research has also been performed on the significance of tool motion parallel to the surface, but these works have continued to assume a constant delay period. As a result the full significance of varying delay lengths is not well understood.

1.3. Research Objective

The goal of this research is to explore some of the effects of the variable delay  period on the stability of a machine tool system in regenerative chatter. Two causes of  variable delay length will be examined: off-axis tool motion and variable speed machining. This will include the lower limit limit of stability and post instability behavior. The problem will be approached from an analytical direction first. A model of the system will be developed that accounts for tool motion along the workpiece and the impact of this motion motion on the delay period. This system will be analyzed to determine its stability.  Next a numerical model of the system will be created that will account more completely for the nonlinearities induced by variation in the delay period. The benefit of  such a model is that it can be used to study the behavior of the system after the onset of 

8 chatter. It can also be used to simulate simulate a particular machining process and compare the result to experimentally collected results. Finally a machine tool fixture will be created that will allow experimental tests to  be conducted for a system whose dynamics are similar to those in the models. The tests will be conducted on a Moore 3 Diamond Turning Machine. This is a very precise and accurate machine tool base with an air bearing spindle. This reduces or eliminates some some sources of experimental error, such as error motion from the spindles or slides. The experimental results will then be compared to the analytical and numerical models of the  process.

1.4. Literature Review

Many different chatter models consider various factors affecting the machining  process. Many of them consider the simplest machining process: orthogonal cutting. An example of this process would be a plunge cut on a lathe using a cut-off tool. While the simple geometric relationships involved make the problem relatively easy to examine analytically this type of process represents only a small portion of turning operations. o perations. True orthogonal cutting is rare in practice as most processes use an oblique geometry in which the tool is at an angle to the cutting velocity. Processes such as milling have more complicated geometric relationships than a plunge cut and are harder to model. However, many processes can be modified into a form that approximates orthogonal cutting and it is therefore a worthwhile starting point for any examination of regenerative chatter.

9 Tlusty [19] presents a variety of machine tool chatter models, one of which models the tool-workpiece system as a spring-mass-damper system as seen in Figure 1.1. In this model only motion motion of the tool in the x-direction is is considered. This is because motion in this direction most directly affects the thickness of the chip and the surface left  behind. For this model the cutting force is modeled as being linearly proportional to the area of the chip being removed. The steady state portion of the cutting force is neglected as it will simply simply result in a static offset. offset. Only the variable part of the cutting force force is considered:  FV 

= Cb( x0 − x)

(1.1) (1.1)

In this equation b is the width of the chip perpendicular to the x the x-direction -direction and the direction of the cutting velocity, C is C is the cutting stiffness or specific cutting energy, which is an empirically derived constant relating the area of the uncut chip to the cutting,  x is the present displacement of the tool from the static load position and x and x is the 0

displacement of the tool on the previous rotation rotation of the workpiece. The spring-massdamper system will have an equation of motion of the form:

&& + c x& + kx = − F sin θ  m x

(1.2)

 Note that as steady state portion of the cutting force is neglected the stability limit is not dependent on the thickness of the chip, only its width. The dependence of the cutting force on the state of the system at a previous time results in what wha t is referred to as a retarded dynamic or delay differential equation. Delay terms can be best dealt with in the frequency domain and most authors approach the problem in this manner. The system can be represented using a block 

10 diagram as seen in Figure 1.2. 1.2. In this figure G(ω) represents the frequency response function of the machine tool structure in the x the x-direction, -direction, F   F  represents the mean cutting m

force, F  force, F  the variable component of the cutting force and exp(j ε) is the delay term, with v

 j = sqrt(-1) and ε being the phase difference between the wave the tool is leaving on the surface and the one it left behind on the previous pass.

m

v

x G( ω   )

-F 

m +

 x

Cb -

exp(j ε)

Figure1.2: Block diagram of cutting forces for regenerative chatter [after Tlusty, 1985].

The block diagram can be analyzed as a feedback feedback control problem. In [19] the  Nyquist criterion is used to determine the point at which the magnitude of the oscillations would neither grow nor decay. This occurs when the open loop transfer function has a value of -1. In this case that is given by: CbG(ω )(1 − e

− jε 

)

= −1

Solving Equation 1.3 for the width of the cut at the stability limit gives:

(1.3)

11 blim

=

−1 CG (1 − e − jε  )

(1.4)

If it is assumed that C is C is purely real then G(1-exp(-jε)) must be real. The minimum value of b of b

lim

is given by: blimmin

=

1

(1.5)

2C Re(Gmin )

It can be shown that for a simple spring-mass-damper system the minimum value of b of b is lim

given by: blimmin

=

2 k  xζ x (1 + ζ x )

(1.6)

C 

This value is the asymptotic value of b of b for stable cutting. This actual limit is higher at certain spindle speeds. speeds. This fact can be seen in a stability chart such as Figure 1.3. This chart was produced for a system with k  = 5500000 N/m, m = 0.48 kg, and ζ = 2.16%. x

x

x

Each of the valleys or ‘lobes’ on the stability chart corresponds to a different number of  chatter undualtions on the machined surface. The first lobe, numbered 1 in Figure 1.3, corresponds to a single oscillation of the tool during each revolution of the workpiece, leaving a single wave on the surface. The next lobe to the left corresponds to two oscillations per revolution and so on. Cutting at width of cut above the lobed line line would  be expected to be unstable and cutting cu tting below should be stable.

12 Stability chart 12 Stability limit, mm 10

 Asy mptotic lim limit it

  m 8   m  ,    t    i   m    i    l 6   y    t    i    l    b    i   a 4    t    S

3

2

2

1

0 0

20000

40000

60000

80000

Spindle Spindle s pee d, RPM RPM

 Figure 1.3: 1.3: Stability chart.

It can be seen from such a diagram that there are certain speeds at which material removal rates can be increased greatly. Unfortunately, for most most machine tools these speeds are above the maximum maximum speed of the machine. Examining the stability chart further it can be seen that the peaks in the stability chart occur when the spindle speed is equal to the system natural frequency (about 538 Hz or 32000 RPM in this system,) divided by an integer, and that the low point of each lobe corresponds to a phase shift of  approximately -90 degrees between waves. These results occur because the frequency at which a system vibrates during regenerative chatter is close to the dominant natural frequency of the system and a phase shift of -90 degrees is the most favorable condition for regeneration. The preceding model accounted for the behavior of a system with a single mode of vibration perpendicular to the workpiece surface. Real systems have multiple modes

13 of vibration which may be oriented in any direction. Tlusty also presents a method method for  dealing with such situations. To deal with the differing orientations orientations of modes of  vibration a directed transfer function is used. These are the projection of the system system transfer function in the x the x-direction. -direction. To account for systems systems with multiple multiple modes of  vibration the stability is based upon the total x total x-direction -direction transfer function of the system which is the sum of the oriented transfer functions of the modes. Ota et. al. [10] also deals with a two degree of freedom system. system. In this case the workpiece is modeled as having two perpendicular modes of vibration, while the tool is considered to be rigid. Although the tool is modeled as having two degrees of freedom, the computation of chip thickness neglects the effect that tool motion parallel to the surface could have on the actual delay period. Experimental measurements of cutting stiffness are shown as are the effects of speed, depth of cut and width of cut on the cutting force. Numerical simulations are also also compared to results from the cutting tests and found to match reasonably well. Some researchers have considered some of the nonlinearities present in machining systems. Nonlinearities can include nonlinear stiffness or damping in the machine tool structure. Another significant nonlinearity is the specific cutting energy. According to Rhorer and Lucca [11] this is not constant, but a function of the chip thickness, with C increasing C increasing as chip thickness decreases. Other machining processes have other types of nonlinearities. For example, in a milling process each flute of the cutter spends only part of its time in contact with the workpiece. To model this a step function must be used.

14 Kato et. al . [6] deals with nonlinear behavior related to the cutting force coefficient. The authors discuss to two types of chatter. The first is a nonregenerative form of chatter referred to to as primary chatter. The second form of chatter is regenerative chatter which is referred to as secondary chatter. This work considers a fixed cutting tool tool cutting a flexible workpiece. This is a reasonable model, as for most machine tools the structure of the tool itself will be very ve ry stiff, but the workpieces may be fairly long and slender. Kato et. al . [7] continues the work begun in [6] by developing a cutting force model that accounts for dynamic effects such as variation in cutting speed, rake angle and cutting depth due to workpiece motion. Experimental results are presented which are claimed to verify the presence of the lag in cutting force described in [6]. Kato, et. al. [8] concludes the research from [6] and [7] by incorporating the  previous results into an analytical model that accounts for the variation of cutting force due to many factors. The workpiece is modeled as a two degree-of-freedom system and the cutting force is affected by many factors including dynamic variations in the chip thickness, cutting velocity and rake angle. The stability of the system is determined and  plotted and the analytical model is compared to experimental results. Hanna and Tobias [3] examined the effect of nonlinear stiffness and nonlinear  specific cutting energy. In this work the stiffness stiffness of the machine tool was modeled as a cubic function of displacement and the specific cutting energy was modeled as a cubic function of chip thickness. It was found that the system system would be conditionally stable. That is to say, if the chatter amplitude exceeded a certain threshold the system would

15  become unstable, but if the system were to vibrate at a lower level the oscillations would decay as seen in Figure 1.4. The system still showed showed an unconditionally stable region, but some of the lobes seen seen in Figure 1.3 were only conditionally stable. This further complicates the process of using the stability lobes to increase material removal rates.

Disturbance amplitude

Unstable

Stable

 blim

Figure 1.4: Conditional stability, [after Hanna and Tobias, 1974].

Another type of nonlinearity is process process damping. This is a phenomenon in which the damping of the system will increase at low cutting speeds. speeds. It is due to plowing and friction generated as the flank face of the cutting tool rubs against the surface of the workpiece at low speeds. The sources of this rubbing are explained well in Elbestawi, et.al. [2]. A model is presented that accounts for for material that is plowed along by the tool rather than being cut off. The end result of process damping damping is that at low cutting speeds the stability limit for for the cutting process will increase. increase. This is demonstrated in Stepan [16]. An additional damping term is added to the structural damping in the system system that is

16 inversely proportional to the spindle speed. Thus damping grows quite large at slow slow speeds. Stepan also considers the impact of a more thorough model of the cutting force. As material is removed from the workpiece surface the cut chip will slide up the rake face of the cutting tool. The friction depends on the length and thickness of the chip sliding up the rake face. Stepan also presents a more thorough thorough model of the cutting force that accounts for the distribution of forces along the rake face of the the tool. This model results in two delay terms being present in the system. The first is the “long delay” which is the length of time time for one rotation of the workpiece. The second is the “short delay” which represents the time required for a chip to slide up the rake face of the tool. Equation 1.7 gives the variable portion of the cutting force as:  F

= C ∫0 w x (ψ )( x( u − d 0π − ψ ) − x( u − ψ )) d ψ    h

(1.7)

In this equation C is C is the conventional cutting force coefficient, h is the length of the chip in contact with the tool, ψis the position along the chip, w x(ψ) is a shape function for the force distribution, and d 0 is the diameter of the workpiece. To incorporate this force model with the rest of the system a change of variable must be performed so that the governing equation of motion is in terms of arc length u rather than time, where u = Vt and V is V is the cutting speed. The result of including both the long and short delays in the model is that the stability of the system may improve slightly at both low and high speeds. Some researchers also examine methods to improve the stability of a machining  process. In addition to making the machine itself stiffer or more highly damped the

17 stability limit for regenerative chatter can be increased by continuously varying the cutting speed. The improvement in stability was explained in Inamura and Sata [5]. This  paper uses an analytical model which uses function space theory and series solutions to study the behavior of the system. The variable rotation period improves stability stability by causing variations in the phase of the present and previous positions of the cutting tool. The authors also predict quite a large increase in the stability limit, in some cases up to a factor of 10. Takemura et. al. [18] presents experimental results from variable speed turning on a lathe. The workpiece was mounted on a flexible arbor so so that the main compliance of  the system was in the workpiece itself rather than the tool. The authors used different wave forms for the speed speed variation including triangular, and rectangular. It was found that if the spindle speed was varied vibration amplitudes could be greatly reduced for  conditions that were unstable during a constant speed cut. The authors suggest that variable speed machining is a practical and useful method to increase stability as it requires no special hardware on most machine machine tools and no sensors. They also state that due to the open loop nature of such a method it will not cause instabilities of its own. This last fact is contradicted by results in Sexton et. al. [14]. This work presents presents a stability model that shows that under some conditions stability can be decreased by the varying speed. These cases seem to mostly mostly involve the upper peak of stability lobes in in constant speed machining which can be clipped off in variable speed machining. machining. Sexton also explained some of the problems with previous analytical studies of variable speed cutting, including [5], which predicted stability gains far larger than were seen in

18 experiments. The weakness of the previous work was that that it was assumed that chatter  only occurred at the natural frequency of the machine tool system and system excitations at other frequencies were neglected. According to Sexton a stability gain of between 50% and 100% can reasonably be obtained. Sexton and Stone [13] present present results of  computer simulations using an analog computer as well as experimental results. Another problem with the use of variable speed machining is explored in Sexton and Stone [12]. During variable speed machining large transient vibrations can occur due to the changing spindle speed even in a stable cut. These vibrations are usually not as severe as they would be in an unstable cut but they make it difficult to obtain good surface finish. As a result variable speed machining is better suited to rough cutting rather than finish cutting. A more recent work on the subject, Tsao et. al. [20] models a system with variable speed using the spindle angular position as the independent variable, instead of  time. The result of this is a time varying system with a constant delay length. This makes it much easier to analyze than a system with fixed coefficients and a variable delay length. The authors claim that the methods they develop are suitable to many other types of machining stability problems. The origin and nature of machining instabilities have been studied by a variety of  authors and a variety of means. Each model considers the importance of different different factors that can impact the stability of the cutting process. Despite the work that has been done there are still many areas that have not been studied adequately and the potential exists to

19 significantly improve machine performance with a more complete understanding of the subject.

1.5. Research Plan

The first phase of this project is the development of a model of orthogonal cutting using a two degree-of-freedom model. This differs from from most existing models models in that the tool is allowed to move along the workpiece in the direction of the cutting velocity, as shown in Figure 1.5. 1.5. This part of the project involves a search of relevant literature to obtain background knowledge about interactions between cutting tools and workpieces, and existing models of chatter for different types of machining processes, such as end milling. A model of chatter in orthogonal orthogonal cutting with two degrees-of-freedom is is then developed using a linear system model and a simple linear model of the cutting force. The models stability is analyzed using the tau-decomposition method and MATLAB  programs were written to produce stability charts.  Next a numerical model of the process is developed. This model simulates the cutting process using an RK45 time domain integration scheme. This program can also create stability charts based upon the simulated results.

20

Cutting force vector 

Chip flow

 N , Spindle speed in RPM

 x

Θ

k  x m x ,  , m y c x c y

k  y

Figure 1.5: Two degree-of-freedom system model.

The next phase of this project consists of the design and construction of a test flexure to replicate the two degree-of-freedom degree-of-freedom model. This was used to perform cutting tests with a Moore 3 Diamond Turning Machine. The flexure has about 10% of the stiffness of the machine tool in the directions normal and tangential to the workpiece surface. The motivation for designing the flexure flexure with deliberately high compliance is so that the critical loop dynamics will be well understood. By making the compliance high, the majority of the motion of the cutter will be due to the simple and well known dynamics of the flexure, rather than the less well known dynamics of the machine tool itself. This also lowers the stability stability limit so that instabilities can be studied on the

21 diamond turning machines used for the experimental experimental testing. The results of cutting tests are compared to the numerical simulation to gauge its accuracy. The numerical simulations is used to study variable speed machining. Experimental results are compared to the simulation. The effect of different spindle spindle oscillation frequencies and amplitudes on the stability of a machining process are examined using the simulation to predict the gain in stability under different conditions. Results are also compared to experimental data from the literature.

22

Chapter 2

ANALYTICAL TWO DEGREE OF FREEDOM MODEL

2.1. Derivation of Equation of Motion

To arrive at a model that can be handled with analytical methods some simplifications were made. First it is assumed assumed that the cutting force is proportional to the area of the chip, and that it acts at a constant angle to the workpiece surface. The stiffness and damping of the machine machine tool are assumed to be linear. As a result the steady deflection caused by the mean cutting force can be neglected. The resulting equations of motion are:  x-direction:  x-direction:

&& + c x x& + k  x x = −α  F  m x x

(2.1)

where: α  = sin Θ

 y-direction:  y-direction: m y &y& + c y y& + k y y = − β F   where:  β  = cos Θ

(2.2)

23 Cutting force approximation:  F  = bC ( x (t ) − x(t − τ ))

(2.3)

Where b is the width of the cut, C is C is a material dependent cutting stiffness term, and τ is τ  is the length of time since the tool tool last passed over this point on the workpiece. For the single degree-of-freedom case this amount is fixed and is given by: τ  = 60 / N 

where N  where N is is the spindle speed in RPM.

2.1.1. Conversion of Equation of Motion to Arc Length Rather Than Time

For the two degree-of-freedom model the magnitude of the time delay is not constant due to yto y-direction direction displacement. Because constant delay terms are much easier  to deal with, the equations of motions are converted to a form in terms of arc length, u, defined as follows: u = Vt Vt + y where V is V is the mean linear speed of the tool relative to the workpiece given by: V

= π dN  / 60

wher wheree d is d is the workpiece diameter 

To convert the equations we need the following: du dt 

= V + y& ≈ V 

This should be a reasonable approximation since V is V is much greater than  y& for typical situations. For example, consider a tool with a y a y-direction -direction natural frequency of 300 Hz, a

24 100 mm workpiece diameter and a spindle speed of 600 RPM. In this situation,V  situation, V will will  be given by Equation 2.4: V

( 60 sec / min) = 3.14 m / se s ec = πdN  / 60 = 0.1 m *π  * 600 RPM / (6

(2.4)

Assuming an amplitude of 10 microns (which would be reasonable during chatter  oscillations) and that the tool moves in a sinusoidal manner  y& will be given by Equation 2.5:  y&

= 2π f 10 microns = 0.0189 m / sec

(2.5)

The equations of motion are written as:  x-direction:  x-direction:

 du  m x    dt  

2

x ' '+ c x

du dt 

x '+ k x x

= −αbC( x( u) − x( u − dπ   ))

(2.6)

 y-direction:  y-direction:

 du  m y    dt  

2

y ' '+ c y

du dt 

y '+ k y y

= − β bC( x( u) − x( u − dπ   ))

(2.7)

Terms such as  x′ ,  y ′ indicate derivatives with respect to arc length rather than time. The delay term is now the constant length d π rather π  rather than a varying length of time.

25 2.1.2. Matrix Form and Characteristic Equation

The equations are then simplified somewhat by dividing the x the x-direction -direction equation through 2

2

2

2

 by m x(du/ du/dt ) = m xV  , and the y the y-direction -direction equation by m y(du/ du/dt ) = m yV  , and defining the following terms: c~ x

=

c~ y

=

c x m xV  c y m yV 

~ k  x

=

~ k  y

=

k  x m xV 2 k  y 2

m yV 

~  F  x

=

~  F  y

=

− α bC  m xV 2

− β bC  m yV 2

Using these new symbols the equation of motion becomes:  x-direction:  x-direction: ~ ~  x ' '+ c~ x x'+ k x x = Fx ( x( u) − x( u − d π ))

(2.8)

 y-direction:  y-direction: ~ ~  y ' '+ c~ y y'+ k y y = Fy ( x( u) − x( u − d π ))

(2.9)

The equations of motion can ca n be written in matrix form as:

1 0 

c~ x  + 1  y' '  0 0  x ' '

~ ~ 0   x'  k x − F x  + ~ c~   y ' − F  y

   

 y

~ 0   x   F  x ~   +  ~ k y   y   F  y

0  x( u − d π ) 

0 =    0  y ( u − d π )  0

(2.10)

Taking the Laplace Transform and determinant to find the characteristic equation yields:  D( λ ) = λ4 ~

~

~

~ 2 + ( c~ y + ~cx )λ3 + ( k y + ~cx ~cy + k x − F  λ x) ~

~ ~

  ~

+ (c~ x k y + c~y ( k x − F~x ))λ + k y ( k x − F~x ) + ( λ2 + ~cy λ  + k y ) F~x e− λd π 

(2.11)

26 2.1.3. -Decomposition Form

For the stability testing method used (presented in [4] and [9],) it is necessary nec essary to rewrite the characteristic equation as follows: Define the following terms: a4

~ = 1 / F   x

a3

~ = ( c~ y + c~x ) / F   x

a1

= ( c~ x k y + ~cy ( k x − F~x )) / F~x

~

~

a2

~

~

~ = ( k y + ~cx c~y + ( k x − F~x ) ) / F   x

 

a0

~ ~

~ = ( k y ( k x − F~x ) ) / F   x

Combining these terms with with (2.11) gives: ~

~  D( λ ) / F x

= a4 λ4 + a3λ3 + a2 λ2 + a1λ + a0 + (λ2 + ~cy λ + k  y ) e− λd π 

(2.12)

Setting D Setting D equal to 0, this becomes: ~

e

λd π 

− ( λ2 + c~y λ  + k y ) = a 4 λ4 + a 3λ 3 + a2 λ 2 + a1λ +  a 0

(2.13)

Equation (2.13) is then separated into two parts: w1

= e λd π 

w2

+ c~ y λ  + k y =− 4 a 4 λ + a 3λ 3 + a2 λ 2 + a1λ +  a 0

(2.14) ~

λ2

(2.15)

Letting λ = λ  = iω  the roots of  D  D will occur when |w |w2(iω )| )| = 1, therefore:

1=

( k~

 y

− ω 2 ) + ic~yω 

(2.16)

( a ω − a ω + a ) + i( a ω − a ω ) 4

2

4

3

2

0

1

3

 

(2.16) is squared resulting in the following:

(a ω − a ω + a ) + ( a ω − a ω ) = ( k~ − ω ) +  ~c ω  4

4

2

2

2

0

3

1

3

2

2

 y

2

2

y

2

(2.17)

27 Expanding this and collecting terms yields: a 42ω 8

+ ( − 2a4 a2 + a32 )ω 6 + (2a4 a0 + a22 − 2 a3 a1 − 1)ω 4 ~

(2.18)

~

+ ( − 2a2 a0 + a12 + 2 k y − ~cy2 )ω 2 + a02 − k y2  = 0 The roots of this equation must must now be found. Since (2.16) was squared to produce (2.18), the number of roots found will be twice the number that are actually present in the system. As a result, only the positive positive real roots of the equation need to be examined. Each positive real root y root y j is substituted back into (2.15) to find w2(iy j). The angle of the resulting number is computed as follows: α  j

= tan −1

ℑ( w2 (iy j ))

(2.19)

ℜ( w2 (iy j ))

The angle is then used to generate a sequence of delay values, τ , where τ  jk

= (α j + 2 kπ ) / y j

(2.20)

where k = k = 0,1,2,3...

One of these sequences is produced for each root. The sequences from all the roots are combined and sorted in increasing order. This sequence represents intervals on the time time delay axis. It is now necessary to determine whether the curve w2 enters or leaves the unit circle at each root y root y j. This is shown in Figure 2.1. This is done by evaluating the magnitude of w of w2 at values slightly larger and smaller smaller than the roots. Every time the curve enters the unit circle the number of positive real roots of the original characteristic equation decreases by two. Every time the curve leaves the circle the the number of roots

28 increases by two. The system is only stable stable in intervals where the number of these roots is zero.

1.25

w1

0.75

0.25

  y   r   a   n    i   g   a   m    I

w2

-0.25

y j

y j -0.75

-1.25 -1.25

-0.75

-0.25

0.25

0.75

1.25

Real

Figure 2.1: Plot of w1, the unit circle and roots in the complex plane.

Each value in the τ sequence τ  sequence corresponds to a root of Equation 2.16 and therefore to a crossing of the unit circle. circle. To determine the stability of the system in each interval of the τ sequence τ  sequence each crossing is examined in order. As each is checked the number of   positive real roots is tracked. Every time the τ value τ  value represents crossing into the circle the number of roots increases by two. Each time the crossing represents the number of  roots decreases by two.

29 In this particular problem, the value of the delay is fixed and known so all that remains is to check the stability stability of the interval in which the the delay lies. This gives the stability of a single point in the parameter space of the problem.

2.2. Stability Chart Generation

To generate more useful information a program was written to find the limiting width of cut b for a given range range of spindle speeds. This was used to produce Figure 2.2. On this plot the straight horizontal line is based on the results of [19], which includes the derivation of the stability limit limit for a single degree-of-freedom degree-of-freedom model. The lobed line was  produced by the two degree-of-freedom model described above. Two degree -of-f -of-freedom reedom stability chart 6 Stability limit 5

 Asy mptotic ptotic limit limit

  m 4   m  ,    t    i   m    i 3    l   y    t    i    l    i    b   a 2    t    S

1 0 0

5000

10000

15000

Spindle Spindle spe ed, RPM RPM

Figure 2.2: Two degree-of-freedom degree -of-freedom stability chart.

20000

30 This method has been used to examine the stability of both one and two degreeof-freedom systems. In every case it was found that the stability of the two degree-offreedom model matches the stability of the single degree-of-freedom system.

2.3. Analytical Model Conclusions

Based upon the results of the analytical model of chatter in a two degree of  freedom system the y the y-direction -direction motion has no effect on the stability limit in orthogonal cutting. However, this model has several limitations. Most important are some of the simplifications needed to arrive at a manageable system model. It is possible that these these could affect the results. In addition this technique was only able to predict the limit of chatter and say nothing about post instability instability behavior. To examine this behavior after chatter a more thorough model that can better deal with the nonlinear nature of the system is required. Such a model would be very difficult to derive and work with analytically. As a result a numerical model of the process was developed.

31

Chapter 3

NUMERICAL CUTTING MODEL

3.1. Characteristics of Numerical Model

To further explore the nature of the two degree-of-freedom cutting model a numerical simulation of the process process was developed. This simulation is able to overcome some of the weaknesses of the analytical model, most importantly the requirement that the system be reduced to a linear approximation approximation of the true system. This also overcomes the restriction that the delay length be constant. The simulation computes computes the chip thickness at each time step based upon the present and previous positions of the cutting tool. By computing chip thickness in this way the program eliminates the need to make assumptions about the length of the delay period. As a result, the simulation is able to predict the time domain response of the system described by Equations 2.1, 2.2 and 2.3 without the need for the simplifications and modifications introduced to obtain an analytical solution. This is the system pictured in Figure 1.5. The numerical integration scheme employed is an adaptive Runge-Kutta 45 method. Runge-Kutta methods compute the state of a system, in this case the velocity and displacement of the cutting tool, at a series of discrete time time steps. The state of the

32 system at any time step is based upon the state of the system at the previous time step, as well as any forcing functions. In this system the forcing function is due to the regenerative effect and as a result the state at each time step can depend upon the position of the tool one or more workpiece rotations in the past. The program computes the next state using both a fourth and fifth order  integration technique. By comparing the two results the program is able to estimate estimate the error in the fifth order integration. Upper and lower limits on the estimated estimated error are set in the program. If the error exceeds the maximum maximum allowable error the time step is reduced and the step is repeated. If the error is less than the minimum allowable error the size of the time step is increased for the subsequent time step. By using an adaptive scheme the program can automatically optimize the size of  each time step. As a result the program can determine the response for a system system using fewer time steps than a fixed step size simulation, speeding up its execution. At each time step the position of the tool relative to the workpiece is stored in memory. If the program is so configured it will write write the time, position, velocity and acceleration of the tool, chip thickness and the position of the tool relative to the workpiece to a file for later use. The program can also be configured to only write write to the data file at certain intervals, for instance every ten time steps. This can greatly reduce the size of the resulting file and increase simulation speed. Storage space for history data is dynamically allocated as the program runs, making it possible to optimize the program to run on computers with different amounts of memory by limiting the number of time steps the program can perform.

33 For ease of use and to help reduce the occurrence of errors in the setup of the  program, a configuration file is read when the program is started. An example of such a file is seen in the appendix. This file lets the the user configure many parameters related to to the execution of the program without the need to recompile the source code each time a change is made. Most importantly the configuration file contains parameters of the system being simulated. These include the stiffness, effective effective mass and damping ratio of  the tool in each direction. Workpiece parameters can also be set in the file, including the diameter, rotational speed and specific cutting energy. The feed rate in meters/revolution is also set in the configuration configuration file. For this particular process geometry this this will equal the average chip thickness. Simulation parameters are configured in the file as well. These include the maximum and minimum allowable error tolerance, the amount by which to increase and decrease the time step by if the error tolerance is not met, and the maximum and minimum time step step size. The total length of the simulation simulation and the maximum number of  time steps can also be set in the configuration file. The program can operate in a number of different modes and the configuration file can be used to set the desired mode. The first mode mode is time domain simulation. In this mode the program simulates the behavior of the cutting tool for a desired length of  time under the conditions specified in the the configuration file. The state of the system at each time step is saved for later use. The program will run until either the desired end time is reached or the program has performed the maximum number of time steps specified in the configuration file.

34 A variety of options can be used when operating in this or any mode. In addition to simulating conventional machining systems the system can simulate systems with continuously varying spindle speeds. Such systems can include conventional lathes in which the feed rate changes with the spindle speed resulting in constant chip thicknesses, and machine tools such as diamond turning machines in which the feed rate is typically independent of the spindle speed. In the second type of system the thickness of the chip being removed will vary with the spindle speed. To simplify the process of analyzing and comparing the results of several simulations under different conditions a mode was created in which the program saves not individual time step data but the peak level for each system state (position, velocity and acceleration) during each revolution of the workpiece. The use of such chatter  envelopes results in much smaller data files that are easier to plot and examine. An example of such a plot and its relationship to the full time history is shown in Figure 3.1. In addition to performing time domain simulations of orthogonal cutting the simulation can produce stability charts by examining the growth or decay of oscillations in the simulated system at a variety of cutting speeds and widths of cut. The program determines whether the oscillations grow or decay for a specified number of workpiece rotations. If they are found to decay a set number of times the system is considered to be stable. If they increase a certain number of times or if the cutting tool leaves the workpiece, the system is considered to be unstable.

35 Comparison of envelope and time domain 0.4    2

Full time domain results

0.3

  s    /   m 0.2  ,   n   o    i    t 0.1   a   r    l   e 0   c   c   a   n -0.1   o    i    t   c   e   r -0.2    i    d        X

Vibration envelope

-0.3 -0.4

0

0.5

1

1.5

2

Time, s

Figure 3.1: Comparison of envelope and full time domain results in an unstable cut. To produce a full stability chart the program begins by examining the stability of  the system at two widths of cut specified in the setup file. file. If the upper limit of this range is found to be unstable and the lower range is found to be stable the program then the limit of stability limit is known to be located within that range. It then checks the stability of the system at the middle of the range. If this is found found to be stable the stability limit lies between the upper end of the original range and the middle. This process of   bisecting the range is repeated until the stability limit has been bracketed within a specified tolerance. This process is then repeated for each spindle speed in the range of  interest. Stability charts can be produced for both constant and variable speed machining machining  processes.

36 3.2. Results of Numerical Simulation

The simulation has been used to predict the behavior of a variety of machine tool systems with a variety of parameters. These systems have included single degree of  freedom systems, systems with identical or similar characteristics in the two directions, and systems with widely differing properties in the x the  x and y and y directions. Some of the cutting processes have been stable and some unstable.

3.2.1. Simulation of Stable Process

The simulation was used to predict the behavior of a process with the following characteristics:  x-direction  x-direction stiffness, MN/m  x-direction  x-direction natural freq., Hz  x-direction  x-direction damping ratio  y-direction  y-direction stiffness, MN/m  y-direction  y-direction natural freq., Hz  y-direction  y-direction damping ratio 2 Workpiece specific cutting stiffness, GN/m Workpiece diameter, m Spindle speed, RPM Depth of cut, m/revolution Width of cut, m Cutting force angle, θ, degrees

5.5 538 2.16% 3.3 332 1.3% 1.0 0.050 600 0.00001 0.0005 20

Table 3.1: Stable process simulation parameters.

Based upon the system properties given here it would be expected that the static tool deflection in the x the x and y and y directions would be:  x-direction:  x-direction:

37

Cbt 0 sin Θ  K  x

=

1 * 10 9 * 0.0005 * 0.00001 * sin( 20°) 5.5 * 10 6

= 3.11 *10 −7 m

(3.5)

 y-direction:  y-direction: Cbt 0 sin Θ  K  y

=

1 * 10 9 * 0.0005 * 0.00001 * cos(20°) 3.3 * 10 6

= 1.42 *10 −6 m

(3.6)

These predictions were confirmed by the results of the simulation plotted in Figures 3 .2 and 3.3. It can be seen that as the tool tool entered the workpiece it deflected the expected amount and the chip thickness quickly approached the specified depth of cut, as seen in Figure 3.4. The oscillations in the tool displacement are too small to be seen compared compared to the steady state offset, but the plot of  x-direction  x-direction acceleration in Figure 3.5 shows that there are some oscillations present and that they quickly decay.

3.2.2. Simulation of Unstable Cutting Process

The damping ratios of the stable system described in the previous section were reduced to 0.9% in the x the x-direction -direction and 0.13% in the y the y-direction. -direction. The resulting cutting  process was unstable. In most unstable systems the amplitude of oscillations will grow until the tool leaves the workpiece resulting in an interrupted cut as is the case with this  particular system. This can be seen in the plots of the system response in Figures Figures 3.6-3.11. Initially the cutting tool is continuously in contact with with the workpiece. At this point the

38 Stable cutting process 0

0.5

1

1.5

2

0   m-0.05  ,    t   n -0.1   e   m   e   c -0.15   a    l   p   s    i    d -0.2   n   o    i    t   c -0.25   e   r    i    d        X -0.3

Workpiece rotation period

-0.35 Time, s

Figure 3.2: Stable cutting process x process x-displacement. -displacement.

Stable cutting process 0

0.5

1

1.5

0 -0.2

  m

 , -0.4    t   n   e   m-0.6   e   c   a    l   p -0.8   s    i    d   n -1   o    i    t   c   e -1.2   r    i    d        Y -1.4

Workpiece rotation period

-1.6 Time, s

Figure 3.3: Stable cutting process y process y-displacement. -displacement.

2

39 Stable cutting process 12 10   m

8

 ,   s   s   e   n    k 6   c    i    h    t   p 4    i    h    C

Workpiece rotation period

2 0 0

0.5

1

1.5

2

Time, s

Figure 3.4: Stable cutting process chip thickness.

Stable cutting process 0.008    2

0.006

  s    /   m 0.004  ,   n   o    i 0.002    t   a   r   e 0    l   e   c   c -0.002   a   n   o -0.004    i    t   c   e   r -0.006    i    d      X

Workpiece rotati on period

-0.008 -0.01 0

0.5

1 Time, s

Figure 3.5: Stable cutting process x process x-direction -direction acceleration.

1.5

2

40 oscillation levels grow with each revolution of the workpiece. Figure 3.12 is the the spectrum of the chip thickness before the cut has become interrupted. Recall that the cutting force is modeled as proportional to to the chip thickness. Therefore the distribution of the cutting force spectrum will be the same as the spectrum of the chip thickness. In Figure 3.12 a single peak can be seen in the spectrum at a frequency slightly above the x the x-direction frequency of the system. The shift in the frequency is due to to the increase in cutting force with increasing chip thickness. This effect results in in an added stiffness term in the system, raising the natural frequency slightly. Roughly 2.5 seconds after the start of the cut the amplitude of the cutting tools motion grows so large that the tool begins to leave the surface of the workpiece, and the chip thickness periodically drops to zero. As the tool is periodically out of the workpiece the surface left by the tool is not a continuous wave form, and the time delay period  between the current position of the tool and the surface it is cutting may span more than one workpiece revolution. In addition there is no force at all on the tool when it is out of  the workpiece and at these times the tool vibrates freely. The result is that the regenerative effect is disrupted, and the amplitude of the oscillations levels out. The spectrum of the chip thickness after the tool has begun leaving the workpiece is seen in Figure 3.13. It can be seen that once the tool begins to to leave the workpiece the spectrum no longer has a single large peak. Instead the spectral content is spread spread out across a wider range of frequencies. In addition the spectrum of the x the x and y and y-direction -direction displacements after the cut has  become interrupted are shown in Figures 3.14 and 3.15. The x The x-direction -direction spectrum

41 Unstable cutting simulation 15

Workpiece rotation period   m 10  ,    t   n   e 5   m   e   c   a    l   p 0   s    i    d   n   o -5    i    t   c   e   r    i    d   - -10    X

-15 0

1

2

3

4

5

6

Time, s

Figure 3.6: Unstable cutting process xprocess x-direction direction displacement.

Unstable cutting simulation 0.4   n 0.2   o   r   c    i 0   m  ,    t   n -0.2   e   m   e   c -0.4   a    l   p   s    i    d -0.6   n   o    i    t   c -0.8   e   r    i    d -1      X

Workpiece rotation period

-1.2 0

0.5

1

1.5

2

Time, s

Figure 3.7: Unstable cutting process xprocess x-direction direction displacement, early growth.

42 Unstable cutting simulation 0.5   m

0

 , -0.5    t   n   e   m -1   e   c   a    l   p -1.5   s    i    d   n -2   o    i    t   c   e -2.5   r    i    d      Y -3

Workpiece rotation period

-3.5 0

1

2

3

4

5

6

Time, s

Figure 3.8: Unstable cutting process yprocess y-direction direction displacement.

Unstable cutting simulation -1.2   n -1.25   o   r   c    i   m -1.3  ,    t   n   e   m-1.35   e   c   a    l   p -1.4   s    i    d   n   o -1.45    i    t   c   e   r    i    d   - -1.5      Y

Workpiece rotation period

-1.55 0

0.5

1

1.5

2

Time, s

Figure 3.9: Unstable cutting process yprocess y-direction direction displacement, early growth.

43 Unstable cutting simulation 30 25

Workpiece rotation period

  m

20

 ,   s   s   e   n    k 15   c    i    h    t   p 10    i    h    C

5 0 0

1

2

3

4

5

6

Time, s

Figure 3.10: Unstable cutting process chip thickness.

Unstable cutting simulation 0.5   m

0

 , -0.5    t   n   e   m -1   e   c   a    l   p -1.5   s    i    d   n -2   o    i    t   c   e -2.5   r    i    d      Y -3

-3.5 -15

-10

-5

0

5

X-direction displacement,

Figure 3.11: Unstable cutting process tool path.

10 m

15

44 Unstable process, cut not yet interrupted 1.E-05 1.E-06   m1.E-07  ,   m   u 1.E-08   r    t   c   e   p 1.E-09   s   s   s 1.E-10   e   n    k   c 1.E-11    i    h    t   p 1.E-12    i    h    C

1.E-13 1.E-14 0

200

400

600

800

1000

Frequency, Hz

Figure 3.12: Unstable cutting process chip thickness spectrum.

Unstable cutting process, steady state 1.E-05 1.E-06

  m  ,   m   u 1.E-07   r    t   c   e   p   s   s 1.E-08   s   e   n    k   c 1.E-09    i    h    t   p    i    h 1.E-10    C

1.E-11 0

200

400

600

800

Frequency, Hz

Figure 3.13: Interrupted cutting process chip thickness spectrum.

1000

45 Unstable cutting process, steady state 1.E-05 1.E-06 1.E-07

  m  ,   n 1.E-08   o    i    t   o 1.E-09   m   n   o 1.E-10    i    t   c   e   r 1.E-11    i    d        X 1.E-12

1.E-13 1.E-14 0

200

400

600

800

1000

Frequency, Hz

Figure 3.14: Interrupted cutting process x process x-direction -direction displacement spectrum.

Unstable cutting process, steady state 1.E-06   m  , 1.E-07   m   u   r    t   c 1.E-08   e   p   s    t   n 1.E-09   e   m   e   c   a 1.E-10    l   p   s    i    d   n 1.E-11   o    i    t   c   e 1.E-12   r    i    d        Y

1.E-13 0

200

400

600

800

1000

Frequency, Hz

Figure 3.15: Interrupted cutting process y process y-direction -direction displacement spectrum.

46

Workpiece surface profile

Workpiece surface Current tool pass Previous tool pass Two tool passes ago

Figure 3.16: Workpiece surface profile, unstable cutting process.

47 contains a single dominant peak near the natural frequency frequency of the system. This is the dominant frequency of the chip thickness and the cutting force and also close to the natural frequency of the xthe x-direction direction mode. The y The y-direction -direction spectrum has two dominant frequencies. The higher frequency matches the x the x-direction -direction frequency and is due to the cutting force. The second lower frequency is the the y  y-direction -direction frequency and is present due to the periods of free vibration when the tool is out of the workpiece. The profile of the workpiece at the end of the cut is shown in Figure 3.16. To make it easier to see the shape, surface variations from the average diameter have been exaggerated by a factor of 400. The workpiece in the figure is rotating clockwise. The  plot also shows portions of the cutting tool path relative to the workpiece on previous workpiece rotations, demonstrating how the surface depends on tool positions during multiple workpiece revolutions.

3.2.3. Simulation of Unstable Non-Interrupted Cutting Process

Oscillations in the unstable system discussed in the previous section grew until the cutting tool began to leave the workpiece. For certain unstable systems it was found that oscillation levels would reach a steady state without the cutting tool leaving the workpiece surface. This is a result of motion motion in the y the y-direction. -direction. This effect was seen primarily in systems with similar  x  x and y and y direction modal frequencies. A practical example of such a system system would be a boring bar used to remove material from the inside inside of a cylinder. The system in this particular particular example had the  properties in Table 3.2.

48 From the plots of the tool response in Figures 3.17-3.20 it can be seen that the amplitude of the tool motion reaches a steady state level without the cutting tool ever  leaving the surface of the workpiece. The spectrum of the chip thickness early in the simulation is shown in Figure Figure 3.21. This spectrum is similar similar to the spectrum in in Figure 3.12 in that it has a single dominant frequency. Figure 3.22 is the spectrum later in the cut, once the vibration levels have stopped growing. It is similar similar to the spectrum from from the interrupted cut in Figure 3.13 in that the spectral content has been spread out across a range of frequencies. In the case of interrupted cutting the regenerative effect was disrupted by the tool leaving the cut. In this case the regenerative effect is is disrupted by the motion of the tool in the y the y-direction -direction which causes variations in the delay period.  x-direction  x-direction stiffness, MN/m  x-direction  x-direction natural freq., Hz  x-direction  x-direction damping ratio  y-direction  y-direction stiffness, MN/m  y-direction  y-direction natural freq., Hz  y-direction  y-direction damping ratio 2 Workpiece specific cutting stiffness, GN/m Workpiece diameter, m Spindle speed, RPM Depth of cut, m/revolution Width of cut, m Cutting force angle, θ, degrees

5.0 503 1.7% 5.0 503 1.7% 1.0 0.035 400 0.0001 0.0007 20

Table 3.2: Unstable noninterrupted process simulation parameters.

49 Unstable noninterrupted cutting process 80   m  ,    t   n   e   m   e   c   a    l   p   s    i    d   n   o    i    t   c   e   r    i    d        X

Workpiece rotation period

60 40 20 0 -20 -40 -60 -80

-100 0

5

10

15

20

Time, s

Figure 3.17: Unstable noninterrupted cutting process xprocess x-displacement. displacement.

Unstable noninterrupted cutting process 200   m

Workpiece rotation period

150 100

 ,    t   n   e 50   m   e   c 0   a    l   p   s -50    i    d   n   o -100    i    t   c   e   r -150    i    d        Y -200

-250 0

5

10

15

Time, s

Figure 3.18: Unstable noninterrupted cutting process yprocess y-displacement. displacement.

20

50 Unstable noninterrupted cutting process 250

200   m  ,   s 150   s   e   n    k   c    i    h    t 100   p    i    h    C

50

Workpiece rotation period 0 0

5

10

15

20

Time, s

Figure 3.19: Unstable noninterrupted cutting process chip thickness.

Unstable noninterrupted cutting process 200   m

150

 , 100    t   n   e   m 50   e   c   a    l 0   p   s    i    d   n -50   o    i    t   c   e -100   r    i    d      Y -150

-200 -80

-60

-40

- 20

0

X-direction displacement,

20 m

Figure 3.20: Unstable noninterrupted cutting process tool path.

40

60

51 Unstable noninterrupted cutting process 1.E-07

  m1.E-08  ,   m   u   r    t   c 1.E-09   e   p   s   s   s   e   n 1.E-10    k   c    i    h    t   p    i    h 1.E-11    C

1.E-12 0

200

400

600

800

1000

Frequency, Hz

Figure 3.21: Unstable noninterrupted cutting process chip thickness spectrum prior to steady state.

Unstable Unstable noninterrupte d cutting proce ss, steady state 1.E-04 1.E-05   m  ,   m1.E-06   u   r    t   c   e   p 1.E-07   s   s   s   e 1.E-08   n    k   c    i    h    t 1.E-09   p    i    h    C 1.E-10

1.E-11 0

200

400

600

800

1000

Frequency, Hz

Figure 3.22: Unstable noninterrupted cutting process chip thickness spectrum at steady state.

52 Unstable Unstable noninterrupte d cutting proce ss, steady state 1.E-04   m1.E-05  ,   m   u   r 1.E-06    t   c   e   p 1.E-07   s    t   n   e 1.E-08   m   e   c 1.E-09   a    l   p   s    i 1.E-10    d   n   o 1.E-11    i    t   c   e   r    i 1.E-12    d        X 1.E-13

0

200

400

600

800

1000

Frequency, Hz

Figure 3.23: Unstable noninterrupted cutting process x process x-direction -direction displacement spectrum at steady state.

Unstable Unstable noninterrupte d cutting proce ss, steady state 1.E-04   m  , 1.E-05   m   u   r    t   c 1.E-06   e   p   s    t 1.E-07   n   e   m1.E-08   e   c   a    l   p 1.E-09   s    i    d   n 1.E-10   o    i    t   c   e 1.E-11   r    i    d        Y 1.E-12

0

200

400

600

800

1000

Frequency, Hz

Figure 3.24: Unstable noninterrupted cutting process y process y-direction -direction displacement spectrum at steady state.

53 To further explore this effect a number of simulations were performed with systems similar to the one above. For Figures 3.25-3.26 a series of systems systems were simulated with a range of different y different y-direction -direction modal modal frequencies. frequencies. All other  characteristics of the system were kept the same. For this particular cutting tool system the length of the wave left on the surface by the x the x-direction -direction motion will be given by Equation 3.3:  N  60 sec/ min  f  x

π d  =

400RPM

π 0.035m = 0.00146m 60 sec/ min* 503 Hz 

(3.7)

The vibration levels in each case reach steady state when the y the y-direction -direction displacement amplitude reaches approximately 0.35 mm. This is roughly 25% of the the length of waves left on the surface by the x-direction oscillation. oscillation. For the systems with the 491 Hz and 539 Hz modal frequencies oscillations quickly grow until the tool leaves the workpiece. The envelope plots only extend until the tool leaves the workpiece for these systems. A similar effect effect can be seen in Figure 3.27-3.28. The systems simulated for for this  plot all have the same modal frequencies, damping ratios and x and x-direction -direction stiffness. However all these systems have different values of  y-direction  y-direction stiffness. The legend indicates the y the y-direction -direction stiffness as a percentage of the x the  x-direction -direction stiffness. stiffness. It can be seen that as the y the y-direction -direction stiffness is increased the steady state amplitude o f motion in the x the x-direction -direction increases until the cut becomes interrupted. interrupted. Once again the amplitudes level out when the y the y-direction -direction displacement amplitude amplitude reaches roughly 0.35 mm. mm. In this case the cutting tool with 200 % stiffness left the workpiece.

54 Envelope plots, unstable noninterrupted cutting, variable y-direction natural frequency 150   m 100  ,    t   n   e 50   m   e   c   a    l   p 0   s    i    d   n   o    i -50    t   c   e   r    i    D   - -100      X

491 Hz 495 Hz 503 Hz 519 Hz 529 Hz 539 Hz

-150 0

5

10

15

20

25

30

Time, s

Figure 3.25: X  3.25: X -direction -direction displacement envelope plots, variable y variable y-direction -direction modal frequency, x frequency, x-direction -direction frequency: 503 Hz.

Envelope plots, unstable noninterrupted cutting, variable y-direction natural frequency 200   m

150 100

 ,    t   n   e 50   m   e   c 0   a    l   p   s    i    d -50   n   o -100    i    t   c   e   r    i -150    D        Y

491 Hz 495 Hz 503 Hz 519 Hz 529 Hz 539 Hz

-200 -250

0

5

10

15

20

25

30

Time, s

Figure 3.26: Y -direction -direction displacement envelope plots, variable y variable y-direction -direction modal frequency, x frequency, x-direction -direction frequency: 503 Hz.

55 Envelope plots, unstable noninterrupted cutting, variable y-direction stiffness, constant frequency. 150   m 100  ,    t   n   e 50   m   e   c   a    l   p 0   s    i    d   n   o -50    i    t   c   e   r    i    D   - -100      X

2% 20% 50% 100% 200%

-150 0

5

10

15

20

25

30

Time, s

Figure 3.27: X  3.27: X -direction -direction displacement envelope plots, variable y variable y-direction -direction stiffness.

Envelope plots, unstable noninterrupted cutting, variable y-direction stiffness, constant frequency. 200   m  ,    t   n   e   m   e   c   a    l   p   s    i    d   n   o    i    t   c   e   r    i    D        Y

0 -200

2% 20%

-400

50% 100%

-600

200%

-800

-1000 0

5

10

15

20

25

30

Time, s

Figure 3.28: Y-direction displacement envelope plots, variable y variable y-direction -direction stiffness.

56 Envelope plots, unstable noninterrupted cutting, variable y-direction y-direction damping ratios. ratios. 150   m 100  ,    t   n   e 50   m   e   c   a    l   p 0   s    i    d   n   o -50    i    t   c   e   r    i    D   - -100      X

0.10% 0.70% 1.30% 1.70% 2.50% 4% 10%

-150 0

5

10

15

20

25

30

Time, s

Figure 3.29: X  3.29: X -direction -direction displacement envelope plots, variable y variable y-direction -direction damping.

Envelope plots, unstable noninterrupted cutting, variable y-direction y-direction damping ratios. ratios. 250 200   m

150

0.10%

 ,    t   n 100   e   m   e 50   c   a    l 0   p   s    i    d   n -50   o    i    t   c -100   e   r    i -150    D        Y

0.70% 1.30% 1.70% 2.50% 4% 10%

-200 -250

0

5

10

15

20

25

30

Time, s

Figure 3.30: Y -direction -direction displacement envelope plots, variable y variable y-direction -direction damping.

57 Another series of simulations were performed with a range o f  y-direction  y-direction damping ratios as shown in Figures Figures 3.29-3.30. It can again be seen that amplitude of   y ydirection motion at steady state is similar for a variety of simulated cutting systems. systems. In this case the tool left the workpiece for the 10% damping ratio case.

3.3. Stability Chart Generation

When used to produce a stability chart the numerical simulation confirms the result of the analytical model that the y the y-direction -direction motion does not affect the stability limit. This is shown in Figure 3.31, which presents stability charts for both a single and a two degree-of-freedom system. For both systems systems the x the x-direction -direction system characteristics

Stability chart, numerical and analytical results 12  Analytical

10

Numerical 2-DOF Numerical 1-DOF

  m 8   m  ,    t    i   m    i    L 6   y    t    i    l    i    b   a 4    t    S

2 0 0

10000

20000

30000

40000

50000

Spindle Speed, Rpm

Figure 3.31: Numerically generated stability charts

60000

70000

80000

58 are the same. It can also be seen that the numerical numerical stability model compares well to the analytical model of Chapter 2. This result can be understood as the y the y-direction -direction motion and x and x-direction -direction motion are  both caused by the same forcing function. As a result there can be no y no y-direction -direction motion unless there is some variation in chip thickness and therefore cutting force due to x to x-direction motion. By this time time however chatter has already begun, so the motion motion cannot halt the onset of chatter.

3.4. Conclusions from Numerical Model

Based upon the results obtained using the numerical simulation of orthogonal cutting y cutting y-direction -direction motion cannot prevent the onset of chatter, but it does have an effect on the growth of chatter and on the characteristics of the system once chatter has begun. In particular the motion can slow the growth of chatter and even limit the maximum amplitude of chatter oscillations.

59

Chapter 4

EXPERIMENTAL WORK 

4.1. Experimental Setup

A two degree-of-freedom flexure was designed and constructed to experimentally verify this model. The stiffness of the flexure is much lower than the stiffness stiffness of the machine tool base. Therefore the dynamic behavior of a tool mounted in the flexure is due mainly to the compliance of the flexure. In this way the behavior of the cutting tool approximates the model used to develop the stability predictions above. The design of the flexure is shown in Figure 4.1. The flexure design uses notch hinges that are explained in Smith and Chetwynd [15]. This type of flexure flexure

Figure 4.1: Final flexure design.

depends upon pairs of closely spaced holes for its compliance. In contrast to

60 other types of flexures which might use beam or leaf elements, the characteristics of  notch flexures are dependent mainly mainly on the spacing between the holes. As most deformation occurs in the thin web of material between the holes manufacturing errors in other parts of the flexure will not have a very large impact on the overall system characteristics. As a result it is is possible to make a notch style flexure using a drill press and bandsaw. Other styles of flexures flexures may need to be produced on a milling milling machine, which is often a more time consuming and expensive process. Another advantage of this type of flexure flexure is its monolithic construction. As it was constructed from a single piece of material there are no joints within the flexure that could potentially introduce nonlinearity. This helps to produce cleaner data from sensors such as accelerometers which can be very sensitive to rattling in joints. The flexure was built and found foun d to have a static stiffness of 5.5 MN/m in the x the x-direction and 3.3 MN/m in the y the y-direction. -direction. To determine dynamic characteristics of the flexure a modal analysis was performed pe rformed using a Kistler 500 N impact hammer, a Kistler  triaxial accelerometer with a range of 50 g and a Hewlett-Packard 35670A Dynamic Signal Analyzer. During the tests the accelerometer was kept fixed and the point of  excitation was moved. This was done to avoid potential problems caused by changing the mass loading of the structure. structure. The set of points used to perform the modal analysis and their locations on the actual part can be seen on plot of the flexure and the undeformed mesh seen in Figures 4.2a and 4.2b. The accelerometer was located at point 6 and oriented to the axis shown in Figure 4.2b. Only the x the x and y and y direction signals were recorded.

61 Undeformed mesh

Side view of flexure 2

1 9

3 5

2

3

4

5

6

10 7

8

4 6

10 7

9 1

8

(a) Mode # 1:332 Hz

(c)

(b) Mode # 2: 537 Hz

(e)

Figure 4.2: Experimentally determined flexure mode shapes

Star Modal analysis software was used to analyze the resulting data. The first mode of vibration corresponds to tool translation in the y the y-direction -direction and has a natural frequency of 332 Hz. This is pictured in Figure 4.2c. The second mode of vibration has

62 a natural frequency of 537 Hz and corresponds to tool translation in the x the x-direction -direction as seen in Figure 4.2d. The damping ratios for these modes were more difficult to determine. Star Modal predicted fairly high damping ratios ratios for the modes: 2.16% and 1.3% in the x the x and y and y direction respectively. Other methods such as a circle fit and the half   power point method predicted damping ratios from 0.5% to 1% or more. The use of the synthesis function on the HP signal analyzer predicted damping ratios of 0.16% and 0.13% in the x the x and y and y directions. To gauge the separation of the x the x and y and y direction modes the accelerometer was  placed on the flexure near the tool mounting location. The impact hammer was then used to excite the flexure in both the x the x and y and y directions. From the results in Figures 4.3 and 4.4 it can be seen that the response of the x the x-direction -direction mode to an x an x-direction -direction excitation is much greater than its response to a y a y-direction -direction excitation. Figures 4.5 and 4.6 show a similar result for the y the y-direction -direction mode. As a result the response of each mode should be largely independent of the other mode. There is a small peak at about 430 Hz in these frequency response functions. The response of this mode is is smaller than the other two modes and this frequency was not observed in the data from cutting tests so this mode does not appear to be significant. It is possible that that this mode is a result result of the clamping system used during the modal testing.

63 Flex ure re sponse ponse at tool mounting, x -direction, -direction, excitation in x-direction 10    N    / 1   g  ,    F    R    F    d 0.1   e   s   a    b   n 0.01   o    i    t   a   r   e    l   e   c 0.001   c    A

0.0001 0

200

400

600

800

1000

Frequency, Hz

Figure 4.3: Flexure frequency response function, response and excitation in x in x-direction. -direction.

Flexure response at tool mounting , response in xdirection, ex citation in y-direction y-direction 1    N    /   g  , 0.1    F    R    F    d   e   s   a 0.01    b   n   o    i    t   a   r   e    l   e 0.001   c   c    A

0.0001 0

200

400

600

800

1000

Frequency, Hz

Figure 4.4: Flexure frequency response function, response in yin y-direction, direction, excitation in  x-direction.  x-direction.

64 Flexure response at tool mounting, response and ex citation in y-direction y-direction 10    N    /   g  , 1    F    R    F    d   e   s   a 0.1    b   n   o    i    t   a   r   e    l   e 0.01   c   c    A

0.001 0

200

400

600

800

1000

Frequency, Hz

Figure 4.5: Flexure frequency response function, response and excitation in y in y-direction. -direction.

Flexure response at tool mounting, response in xdirection, ex citation in y-direction y-direction 1    N    / 0.1   g  ,    F    R    F    d 0.01   e   s   a    b   n 0.001   o    i    t   a   r   e    l   e   c 0.0001   c    A

0.00001 0

200

400

600

800

1000

Frequency, Hz

Figure 4.6: Flexure frequency response function, response in x in x-direction, -direction, excitation in  y-direction.  y-direction.

65

4.2. Experimental Results

A high speed steel tool mounted in this flexure was used to make a series of  orthogonal cuts using a Moore 3 jig borer base with an Excello dual-hemisphere air   bearing spindle. The workpiece used was an aluminum shaft 35 mm in diameter that had  been machined on a lathe to create a series of disks protruding from the surface  perpendicular to the axis of rotation as seen in Figure 4.7. The disks were of varying widths so that a range of processes could be explored.

Figure 4.7: Example workpiece design.

To measure the response of the flexure a Kistler triaxial accelerometer was used along with an ADE capacitance probe. An HP signal analyzer was used to perform time captures during the cut. These data were later analyzed and plotted.

66 4.2.1. Stable Cutting Test Results

This cut involved a disk 0.5 mm mm wide at a spindle speed speed of 199.4 RPM. The feed -6

rate, which determines the chip thickness, was 3.47*10 m/revolution. Should the measured x measured x-direction -direction displacement ever exceed the feed per revolution in the negative x negative x-direction the tool will leave the workpiece. As seen in Figure 4.8 the system system quickly reaches a steady state with low vibration levels, and the tool remains in contact with the workpiece at all times. Although the theory suggests that the vibration vibration levels should drop to zero this was not found to be the case. It is possible that the small vibrations are due to spindle runout, slideway error, nonuniformities in the workpiece such as inclusions or  voids, variations in the thickness of the disk being cut, or imbalance in the rotating parts of the system. Although the amplitude of the oscillations is very small, the frequency frequency at which they occur is high. As a result the accelerations associated with these displacements is fairly high as seen in Figures 4.10 and 4.12. The simulation was then configured to replicate replicate this cutting process. process. The damping ratios used were 0.9% in the x the x-direction -direction and 0.4% in the y the y-direction. -direction. It can be seen from Figure 4.9 that the displacement predicted by the simulation matches with the experimental results. Since the simulation does not include the sources of random noise  present in experimental testing the vibrations levels do drop to zero. As a result, the acceleration levels are also very small, much smaller than seen in the experimental testing as seen in Figures 4.11 and 4.13.

67 Stable cut, cut, e xperimental data 0.05 0

  m

 ,    t   n -0.05   e   m   e   c -0.1   a    l   p   s    i    d -0.15   n   o    i    t   c -0.2   e   r    i    d        X -0.25

-0.3 0

2

4

6

8

10

12

14

Time, s

Figure 4.8: X  4.8: X -direction -direction displacement, experimental data, stable cut.

Stable cut, simula simula tion 0.05   m

0

 ,    t   n -0.05   e   m   e   c -0.1   a    l   p   s    i    d -0.15   n   o    i    t   c -0.2   e   r    i    d        X -0.25

-0.3 0

2

4

6

8

10

Time, s

Figure 4.9: X  4.9: X -direction -direction displacement, simulation, stable cut.

12

14

68

Stable cut, cut, e xperimental data 20    2

15

  s    /   m 10  ,   n   o    i    t   a 5   r   e    l   e 0   c   c   a   n -5   o    i    t   c   e -10   r    i    d        X -15

-20 0

2

4

6

8

10

12

14

Time, s

Figure 4.10: X  4.10: X -direction -direction acceleration, experimental data, stable cut.

Stable cut, cut, simulation 6.E-04    2

4.E-04

  s    /   m  ,   n 2.E-04   o    i    t   a   r   e 0.E+00    l   e   c   c   a -2.E-04   n   o    i    t   c -4.E-04   e   r    i    d        X -6.E-04

-8.E-04 0

2

4

6

8

10

Time, s

Figure 4.11: X  4.11: X -direction -direction acceleration, simulation, stable cut.

12

14

69 Stable cut, cut, e xperimental data 8    2

6

  s    /   m 4  ,   n   o    i    t 2   a   r   e    l   e 0   c   c   a   n -2   o    i    t   c   e -4   r    i    D        Y -6

-8 0

2

4

6

8

10

12

14

Time, s

Figure 4.12: Y -direction -direction acceleration, experimental data, stable cut.

Stable cut, cut, simulation 3.E-03    2

2.E-03

  s    /   m  ,   n 1.E-03   o    i    t   a   r   e 0.E+00    l   e   c   c   a -1.E-03   n   o    i    t   c -2.E-03   e   r    i    d        Y -3.E-03

-4.E-03 0

2

4

6

8

10

Time, s

Figure 4.13: Y -direction -direction acceleration, simulation, stable cut.

12

14

70

4.2.2. Pulsed Unstable Cutting Test Results

The process parameters for this cut were similar to the parameters above with the -6

exceptions that the spindle speed is now 600 RPM and the feed rate is now 1.15*10

m/revolution. This cutting process was found to be unstable. From the experimental results of Figures 4.14, 4.16, and 4.18 it can be seen that the vibrations grow until the tool  begins to leave the workpiece. Furthermore the signal appears to pulse every three to four seconds with vibration levels dropping to to nearly zero before rising again. The x The x-direction acceleration and displacement spectra from the experimental data in Figures 4.20 and 4.24 both have a single dominant frequency slightly above the flexure frequency in the x the x-direction. -direction. The y The y-direction -direction acceleration spectrum in Figure 4.22 has two peaks, one near the x the x-direction -direction natural frequency and one near the y the y-direction -direction natural frequency. It is believed that the lower of these two peaks is due to the free vibration period when the tool is not engaged in the workpiece. A simulation of this system produced a similar pulsing vibration signal, with the vibration levels rising and falling approximately every five seconds as seen in Figures 4.15, 4.17 and 4.19. The plot of chip thickness in Figure 4.25 shows that during the  periods when the vibrations drop near zero the tool is out of the workpiece for a roughly half a second. During the periods of large amplitude vibration the cutting tool vibrates so much that it removes a considerable thickness thickness of the workpiece. This reduces the chip thickness and cutting force on following passes. As a result tool vibration vibration levels drop until the tool no longer cuts the workpiece at all. Since the cutting force is zero at this

71  point vibration levels drop further. The tool does not contact the workpiece aga in until the feed moves the tool closer to to the surface. This is seen in Figure 4.26 which plots plots the surface profile, and the position and path of the cutting tool at the beginning of the first such gap. Displacements and surface waviness have both been greatly greatly exaggerated in this plot. Figure 4.21 shows the spectrum of the x the x-direction -direction acceleration based upon the simulation. This is similar to to the spectrum from the experimental data. Similarly in Figure 4.23 the acceleration spectrum in the y the y-direction -direction is comparable to the experimental spectrum of Figure 4.22. The damping ratio used for this simulation was 0.175% in the x the  x-direction -direction and 0.13% in the y the y-direction. -direction. These values are very close to the values determined with the signal analyzer. Although these values are lower than those used to simulate the stable system in the previous example, it is possible that there was more damping in the system in that case due to process damping which would be more likely to occur at low speed. Very small changes in the parameters used in the simulation (such as modal frequency, stiffness and damping) can produce dramatic changes in the rate and amplitude of pulsing. For example, changing the x the x-direction -direction modal frequency to 536.3 Hz and the damping to 0.22% doubles the rate of pulsing. It is likely that small errors in the measurement of system parameters are responsible for the differences between simulated and experimental results.

72 Pulsed vibration signal, experimental data 40   m

30

 , 20    t   n   e   m 10   e   c   a    l   p 0   s    i    d   n -10   o    i    t   c   e -20   r    i    d        X -30

-40 0

2

4

6

8

10

12

14

Time, s

Figure 4.14: X  4.14: X -direction -direction displacement, experimental data, pulsed unstable cut.

Pulsed Pulsed vibration signal, signal, simula tion 40   m

30

 , 20    t   n   e   m 10   c   e   a    l   p 0   s    i    d   n -10   o    i    t   c   e -20   r    i    d        X -30

-40 0

2

4

6

8

10

12

Time, s

Figure 4.15: X  4.15: X -direction -direction displacement, simulation, pulsed unstable cut.

14

73 Pulsed vibration signal, experimental data 400    2

  s 300    /   m  ,    l   e 200   v   e    l   n 100   o    i    t   a   r 0   e    l   e   c   c -100   a   n   o    i    t -200   c   e   r    i    d   - -300      X

-400 0

2

4

6

8

10

12

14

Time, s

Figure 4.16: X  4.16: X -direction -direction acceleration, experimental data, pulsed unstable cut.

Pulsed Pulsed vibration signal, signal, simula tion 400    2

300

  s    /   m 200  ,   n   o    i    t   a 100   r   e    l   e 0   c   c   a   n -100   o    i    t   c   e   r    i -200    d        X -300

-400 0

2

4

6

8

10

12

Time, s

Figure 4.17: X  4.17: X -direction -direction acceleration, simulation, pulsed unstable cut.

14

74 Pulsed vibration signal, experimental data 8    2

6

  s    /   m 4  ,   n   o    i    t   a 2   r   e    l   e 0   c   c   a   n -2   o    i    t   c   e -4   r    i    d        Y -6

-8 0

2

4

6

8

10

12

14

Time, s

Figure 4.18: Y -direction -direction acceleration, experimental data, pulsed unstable cut.

Pulsed Pulsed vibration signal, signal, simula tion 0.4   m

0.2

0  ,    t   n   e   m-0.2   e   c   a    l   p -0.4   s    i    d   n -0.6   o    i    t   c   e -0.8   r    i    d        Y -1 -1.2 0

2

4

6

8

10

12 12

Time, s

Figure 4.19: Y -direction -direction acceleration, simulation, pulsed unstable cut.

14

75 Pulsed vibration signal, experimental data    2

1.E+02

  s    /   m1.E+01  ,   m   u 1.E+00   r    t   c   e   p 1.E-01   s   n   o    i 1.E-02    t   a   r   e    l   e 1.E-03   c   c   a   n 1.E-04   o    i    t   c   e 1.E-05   r    i    d        X 1.E-06

0

200

400

600

800

1000

Frequency, Hz

Figure 4.20: X  4.20: X -direction -direction acceleration spectrum, experimental data, pulsed unstable cut.

Pulsed Pulsed vibration signal, signal, simula ted da ta    2

1.E+02

  s    /   m1.E+01  ,   m   u 1.E+00   r    t   c   e   p 1.E-01   s   n   o    i 1.E-02    t   a   r   e    l   e 1.E-03   c   c   a   n 1.E-04   o    i    t   c   e 1.E-05   r    i    d        X 1.E-06

0

200

400

600

800

1000

Frequency, Hz

Figure 4.21: X  4.21: X -direction -direction acceleration spectrum, simulation, pulsed unstable cut.

76 Pulsed vibration signal, experimental data    2

1.E+01

  s    /   m  , 1.E+00   m   u   r    t   c 1.E-01   e   p   s   n   o    i 1.E-02    t   a   r   e    l   e   c   c 1.E-03   a   n   o    i    t   c 1.E-04   e   r    i    d        Y

1.E-05 0

200

400

600

800

1000

Frequency, Hz

Figure 4.22: Y -direction -direction acceleration spectrum, experimental data, pulsed unstable cut.

Pulsed Pulsed vibration signal, signal, simula ted da ta 1.E+00  ,   m   u 1.E-01   r    t   c   e   p   s   n 1.E-02   o    i    t    2   a   s   r    /   e    l   m   e   c 1.E-03   c   a   n   o    i    t   c   e 1.E-04   r    i    d        Y

1.E-05 0

200

400

600

800

1000

Frequency, Hz

Figure 4.23: Y -direction -direction acceleration spectrum, simulation, pulsed unstable cut.

77 Pulsed vibration signal, experimental data 1.E-05 1.E-06   m  ,   m   u   r 1.E-07    t   c   e   p   s1.E-08    t   n   e   m 1.E-09   e   c    l   p   s1.E-10    i    d   n   o1.E-11    i    t   c   e   r    i 1.E-12    d        X 1.E-13 0

200

400

600

800

1000

Frequency, Hz

Figure 4.24: X  4.24: X -direction -direction displacement spectrum, experimental data, pulsed unstable cut.

Pulsed Pulsed vibration signal, signal, simula tion 12 10   m  , 8   s   s   e   n 6    k   c    i    h    t   p 4    i    h    C

2 0 0

2

4

6

8

10

Time, s

Figure 4.25: Chip thickness, simulation, pulsed unstable cut.

12

14

78

Workpiece surface profile Workpiece surface Tool path

Figure 4.26: Workpiece surface and tool path at end of vibration pulse.

79

4.2.3. Unstable Cutting Test Results

Although in the above cases the simulated and experimental results were found to agree fairly well this was not always the case. A cutting test was performed performed with  parameters similar to those above with the exceptions that the spindle speed was 405 -6

RPM and the feed rate was now 1.73*10 m/revolution. The experimental data in Figures 4.27, 4.29, 4.31 show that the oscillation levels grow quickly and remain at a reasonably constant level throughout the test period. This system was simulated using damping ratios of 0.25% in the x the  x direction and 0.13% in the y the y direction and the results results are plotted in Figures 4.28, 4.30, 4.32. Although the maximum oscillation levels reached are a re similar to those seen in the experimental tests the oscillations seem to grow more slowly and an d at the end the amplitude seems to pulse in a manner similar to that that seen in the previous example. It is possible that the discrepancy is due to some unmodeled factor in the process, such as process damping, or nonlinear  specific cutting energy that becomes more significant significant at lower cutting speeds. It may also be the case that the experimental test was not long enough for the system to reach steady state, and that with a longer test time the system would behave more like the simulation.

80 Unstable cutting process, experimental data 25   m

20 15

 ,    t   n 10   e   m   e 5   c   a    l   p 0   s    i    d   n -5   o    i    t -10   c   e   r    i -15    d        X

-20 -25

0

2

4

6

8

10

12

14

Time, s

Figure 4.27: X  4.27: X -direction -direction displacement, experimental data, unstable cut.

Unstable cut, simulation 25 20   m

15

 ,    t   n 10   e   m   e 5   c   a    l   p 0   s    i    d   n -5   o    i    t -10   c   e   r    i    d   - -15      X

-20 -25

0

2

4

6

8

10

Time, s

Figure 4.28: X  4.28: X -direction -direction displacement, simulation, unstable cut.

12

14

81

Unstable cutting process, experimental data 400 300

   2

  s    /   m 200  ,   n   o    i    t   a 100   r   e    l   e 0   c   c   a   n -100   o    i    t   c   e -200   r    i    d        X -300

-400 0

2

4

6

8

10

12

14

Time, s

Figure 4.29: X  4.29: X -direction -direction acceleration, experimental data, unstable cut.

Unstable cut, simulation 400    2

300

  s    /   m 200  ,   n   o    i    t   a 100   r   e    l   e 0   c   c   a   n -100   o    i    t   c   e -200   r    i    d        X -300

-400 0

2

4

6

8

10

Time, s

Figure 4.30: X  4.30: X -direction -direction acceleration, simulation, unstable cut.

12

14

82 Unstable cutting process, experimental data

6.5    2

  s 4.5    /   m  ,   n 2.5   o    i    t   a   r   e 0.5    l   e   c   c   a   n -1.5   o    i    t   c -3.5   e   r    i    d        Y -5.5

-7.5 0

2

4

6

8

10

12

14

Time, s

Figure 4.31: Y -direction -direction acceleration, experimental data, unstable cut.

Unstable cut, simulation 8    2

6

  s    /   m  ,   n   o    i    t   a   r   e    l   e   c   c   a   n   o    i    t   c   e   r    i    d        Y

4 2 0 -2 -4 -6 -8

-10 -12 0

2

4

6

8

10

Time, s

Figure 4.32: Y -direction -direction acceleration, simulation, unstable cut.

12

14

83

4.3. Conclusions of Experimental Work 

The comparison of the experimental and numerical results provides confidence in the accuracy of the numerical simulation. That it is able to capture some some of the complex  behavior of the pulsing chatter signal implies that it is accurately modeling much of the significant detail in the system. Improvements in the accuracy of the simulation could be obtained by refining the simulation mode. One possible refinement would be a more more detailed model of the cutting  process itself including any variation of specific cutting energy with chip thickness or  cutting speed. Another area where more more detailed modeling could be beneficial would be to account for the effect e ffect of process damping in the system model.

84

Chapter 5

VARIABLE SPEED MACHINING

5.1. Background

In the preceding chapters it was shown that tool motion in the y the y-direction -direction does not affect the stability of a cutting process p rocess as the ythe y-direction direction motion amplitude will only  become significant once the process has become unstable. However, stability can be improved by intentionally introducing a similar effect into a system by continuously varying spindle speed during a cut. This has two effects. The first is to to continuously alter the length of the time time delay  period, disrupting the regenerative effect. The second is to change the speed at which surface waviness left by previous passes moves past the tool. Since the cutting tool oscillates at a fairly constant frequency during a cut, the wavelength of the marks left on the surface will depend on the speed at which the the surface moves past the cutting tool. If  the spindle speed is constant the wavelength will be constant and the frequency at which the tool passes through those waves on subsequent passes will be the same as the frequency at which they were created. If the spindle speed is not the same from one pass to the next the cutting tool will pass through those waves at a different frequency than the tool natural frequency further reducing their ability to excite the tool.

85 5.2. Approach

Most previous studies of variable speed cutting have involved machine tool systems in which the same motor turns the spindle and moves the slide carrying the cutting tool. As a result when the spindle speed is changed the feed rate will also change and the thickness of the chip being removed from the workpiece will remain constant. Modern machine tools, including diamond turning machines, use separate drive systems to control the spindle spindle and move the cutting tool. On most machines of this type the controller can vary the feed rate along with the spindle speed to maintain constant chip thickness. However, this feature is not always present present or used. As a result the feed and therefore the thickness of the chip being removed can vary with the spindle speed. To determine the importance of this variation in chip thickness on stability the simulation program introduced in Chapter 3 includes the capability to simulate and determine the stability of variable speed machining systems, both those with constant and variable feed rates. The stability determination procedure is slightly different from that used for constant speed systems. Rather than looking for peak oscillation levels every rotation of the workpiece, it determines the peak oscillation levels during one speed oscillation cycle.

5.3. Variable Speed Machining Results

The unstable cutting process of Chapter 3.2.2 was simulated simulated again. This time the spindle speed was varied in a sinusoidal manner at a frequency of 0.25 Hz, with an

86 amplitude that was 20% of the mean mean spindle speed. In this case that means that the spindle speed varied from 480 RPM to 720 RPM. The first simulation performed varied the feed rate a long with the spindle speed. As a result the chip thickness was constant. From Figures 5.1-5.3 it can be seen that that the variation in spindle speed has stabilized the previously unstable cutting process. The system was then simulated simulated with a constant feed rate. As seen in Figures 5.45.6, the high frequency oscillations present in the original constant speed example are gone and the oscillation levels levels are relatively static. static. Because the spindle speed is is constantly varying while the feed rate remains constant the feed per workpiece revolution is constantly changing as well. As a result the chip thickness thickness varies with spindle speed, causing variations in the cutting force. This results in some some tool motion, but the amplitude of these oscillations is very small compared to the constant speed system in which the cut actually became interrupted. As diamond turning machines are used to produce parts with extremely strict surface finish requirements and tight tolerances even these small oscillation levels would likely be unacceptable on finished parts. It is still possible that variable spindle speed machining could be used to increase material removal rates during roughing cuts earlier  in the production process.

87 Variable speed ma chining, chining, variable feed rate 0   m-0.05  ,    t   n -0.1   e   m   e -0.15   c   a    l   p   s    i    d -0.2   n   o    i    t   c -0.25   e   r    i    d     x -0.3

-0.35 0

2

4

6

8

10

Time, s

Figure 5.1: X  5.1: X -direction -direction displacement, variable spindle speed, variable feed rate (constant feed/revolution.)

Variable speed ma chining, chining, variable feed rate 0 -0.2

  m

 , -0.4    t   n   e   m-0.6   e   c   a    l   p -0.8   s    i    d   n -1   o    i    t   c   e -1.2   r    i    d     y -1.4

-1.6 0

2

4

6

8

10

Time, s

Figure 5.2: Y -direction -direction displacement, variable spindle speed, variable feed rate (constant feed/revolution.)

88 Variable speed ma chining, chining, variable feed rate 12 10   m

8

 ,   s   s   e   n    k 6   c    i    h    t   p 4    i    h    C

2 0 0

2

4

6

8

10

Time, s

Figure 5.3: Chip thickness, variable spindle speed, variable feed rate (constant feed/revolution.)

Variable speed machining, constant feed rate 0 -0.05

  m

 , -0.1    t   n   e   m-0.15   e   c   a    l   p -0.2   s    i    d   n -0.25   o    i    t   c   e -0.3   r    i    d     x -0.35

-0.4 0

2

4

6

8

10

Time, s

Figure 5.4: X  5.4: X -direction -direction displacement, variable spindle speed, constant feed rate (variable feed/revolution.)

89 Variable speed machining, constant feed rate 0 -0.2

  m

-0.4

 ,    t   n   e -0.6   m   e   c -0.8   a    l   p   s -1    i    d   n   o -1.2    i    t   c   e   r -1.4    i    d     y

-1.6 -1.8 0

2

4

6

8

10

Time, s

Figure 5.5: Y -direction -direction displacement, variable spindle speed, constant feed rate (variable feed/revolution.)

Variable speed machining, constant feed rate 14 12   m10  ,   s   s 8   e   n    k   c    i 6    h    t   p    i    h 4    C

2 0 0

2

4

6

8

10

Time, s

Figure 5.6: Chip thickness, variable spindle speed, constant feed rate (variable feed/revolution.)

90 The simulation was used to generate gene rate stability charts for this system using both constant and variable feed rates. From Figure 5.7 it can be seen that there is no difference in the stability limit limit for constant and variable feed rate systems. The variation of spindle speed increased the stability of the system over the constant speed case at lower spindle speeds. At higher speeds however, the variable speed can alter favorable  phase relationships that result in the high peaks on the stability lobes, reducing the stability at some speeds.

Comparison Comparison of constant constant and variable feed 0.0035 Variable feed Constant f eed

0.003

Constant speed 0.0025   m  ,    t    i 0.002   m    i    l   y    t    i    l    i 0.0015    b   a    t    S

0.001

0.0005 0 0

1000

2000

3000

4000

5000

6000

Spindle Spindle s pee d, RPM RPM

Figure 5.7: Comparison of constant and variable speed cutting stability.

The effect of variable speeds was further explored by running cases with different speed variation frequencies and amplitudes. amplitudes. The properties of the single degree-offreedom system used are shown in Table 5.1.

91 X-direction stiffness, MN/m X-direction natural freq., Hz X-direction damping ratio 2 Workpiece specific cutting stiffness, GN/m Workpiece diameter, m Depth of cut, m/revolution Cutting force angle, θ, degrees

3.0 180 1.3% 1.0 0.020 0.0001 20

Table 5.1: Stable process simulation parameters. First the frequency of the spindle speed variation was held constant over a range of amplitudes. Next the amplitude of the variation variation was kept constant over a range of  frequencies. From the results of Figures 5.8 and 5.9 it can be seen that the gain in stability increases with with the amplitude and frequency of spindle speed variation. For any real machine tool system there will be practical limits on the rate and amount that the spindle speed can be varied, placing some of these higher gains out of reach. The variable speed stability charts are somewhat ‘lumpy,’ particularly at higher  spindle speeds. This results from the stability testing testing method used. If the program were set up to more closely examine the stability at each point this effect would be reduced. Unfortunately, this would substantially increase program execution time. The stability chart in Figure 5.10 was generated for comparison to experimental results from Sexton [13], presented in Table 5.2.

92

Spindle speed, RPM 172 234 322 421 561 749

Stability gain, %, from Sexton [13] 35 46 36 42 75 18

Table 5.2: Experimental results from Sexton [13]

The system used in this has a 50 Hz natural frequency with a 6% damping ratio, and spindle speed is varied at 0.25 Hz with and amplitude of 20% of the mean spindle speed. Sexton determined the the gain in stability at several spindle speeds. The simulation was configured to reproduce reproduce this situation. The gain in stability stability over the constant speed case is presented in Figure 5.11. It can be seen that the experimentally determined gains in stability as reported by Sexton are close to the values predicted by the simulation at most data points. The correspondence is worse at lower speeds which could be due to the impact of process damping which is not included in the simulation. The results also differ significantly at the 561 RPM data point. The reason for this is not known. It is  possible that the data is an outlier in the experimental data, or that there is a problem with the stability analysis of the situation. An experimental cutting test was performed using a single degree-of-freedom flexure similar to the the one presented in Chapter 4. The flexure had a 170 Hz natural frequency, with a stiffness of 3.5 MN/m and a damping ratio ratio of 0.66 %. The tool was

93 Variable Spindle Speed Stability Chart, 0.25 Hz Constant Speed

0.002

5% Speed Variation

0.0018

10% Speed Variation

0.0016

15% Speed Variation 0.0014

  m  ,    t    i 0.0012   m    i    L 0.001   y    t    i    l    i 0.0008    b   a    t    S 0.0006

20% Speed Variation 25% Speed Variation 30% Speed Variation

0.0004 0.0002 0 0

500

1000

1500

2000

Me an Spindle Spindle Speed, RPM RPM

Figure 5.8: Variable speed stability limit comparison, constant frequency.

Variable Spindle Speed Stability Chart, 20 % Speed Variation 0.0018

Constant Speed

0.0016

0.05 Hz

0.0014

0.10 Hz 0.20 Hz

  m0.0012  ,    t    i   m 0.001    i    L   y 0.0008    t    i    l    i    b   a 0.0006    t    S

0.30 Hz 0.40 Hz

0.0004 0.0002 0 0

500

1000

1500

2000

Me an Spindle Spindle Spee d, RPM RPM

Figure 5.9: Variable speed stability limit comparison, constant amplitude.

94 Numerically Predicted Stability Charts, Constant and Variable Speed 1.2E-03

Constant Speed Variable Speed

1.0E-03   m 8.0E-04  ,    t    i   m    i    L 6.0E-04   y    t    i    l    i    b   a 4.0E-04    t    S

2.0E-04 0.0E+00 0

200

400

600

800

Spindle Speed, RPM

Figure 5.10: Variable speed stability limit.

Gain in Stability, Numerical vs Experimental Results Results From Sexton 120

Stability Stability Gain, %

100

Gain from Sexton

80

   %  ,   y    t    i 60    l    i    b   a    t    S 40   n    i   n 20    i   a    G 0

-20 -40 0

100

200

300

400

500

600

700

800

Spindle Speed, RPM

Figure 5.11: Variable speed stability gain, simulated and experimental from Sexton [13].

95 used to perform an orthogonal cut with a width of cut of 0.635 mm. mm. As the start of the cut the spindle speed was held constant at 500 RPM. The tool chattered violently, and the chip removed from the workpiece was discontinuous. Starting approximately 9 seconds later the spindle speed was varied between 390 and 600 RPM every seven seconds. From the plot of the the results in Figure 5.12 it can be seen that the oscillation oscillation amplitude dropped at this point. Although there were still still oscillations present they were not as severe as they had been. Figure 5.13 shows the expected stability chart for constant and variable speed machining under these conditions. From this it would be expected that the system would be stable with the variable spindle speed. From this it is is apparent that the system used is is not adequately modeled by the simulation, either because the simulation is flawed or because the test rig used does not Experimental data, variable speed cutting 150

  s   n   o   r   c    i   m  ,    t   n   e   m   e   c   a    l   p   s    i    d    l   o   o    T

100 50 0 -50 -100 -150 0

5

10

15

Time, s

Figure 5.12: Tool displacement, variable spindle speed after 9 seconds.

20

96 Varia ble speed stabili stabili ty chart chart 2 1.8

Variable speed

1.6

Constant speed

1.4

  m   m  , 1.2    t    i   m    i    l 1   y    t    i    l    i 0.8    b   a    t    S 0.6

0.4 0.2 0 0

200

400

600

800

1000

Spindle Spindle spe ed, RPM RPM

Figure 5.13: Constant and variable spindle speed stability.

 behave as expected. Unfortunately the machine tool available for these tests is not well suited to further testing of variable spindle speed machining.

5.4. Variable Speed Machining Conclusions

Based upon the results of the simulation it is believed that variable speed machining can be used successfully to increase the stability of machining processes both on conventional machine tools like lathes and on tools like diamond turning machines. The benefits of speed variation are greatest at lower spindle speeds, and the variation can  be detrimental at higher spindle speeds by reducing the peaks of the stability lobes.

97

Chapter 6

CONCLUSIONS

6.1. Conclusions

A new two degree-of-freedom model of orthogonal cutting including time delay length variation was developed and studied in a variety of ways. The stability of a cutting process was examined using both analytical and numerical means and the results were found to match. The predictions of the numerical simulation simulation were also compared to experimental data collected using a flexure designed to behave like the system in the model. The numerical simulation was found to predict the behavior seen experimentally. It was shown by both analytical and numerical means that the inclusion of the time delay variation has no significant effect on the critical width of cut for the onset of  chatter. It does however change the behavior of the system after the onset of chatter and can limit the amplitude of chatter cha tter oscillations. The use of variable spindle speed machining to increase stability in cutting  processes was explored using a numerical simulation of the process. It was found that this technique could be used to effectively increase stability in certain situations. Results

98 obtained with the simulation were found to be consistent with experimental results obtained from the literature.

6.2. Future Work 

Future work on the two degree of freedom model would include extending it to other machining systems such as milling. In such systems the tool is constantly oscillating due to the nature of the process so the constant delay assumption could lead to  poor predictions of stability. Experimental verification of the unstable noninterrupted cutting situation could be obtained by performing cutting tests with a boring bar or other o ther similar tool with two closely matched natural frequencies. It would also be useful to experimentally determine the stability limit for variable spindle speed machining. Although this technique has been used in manufacturing situations already, few experimental results have been published. Doing this would make it possible to further gauge the accuracy of the simulation.

99

Appendix A

SIMULATION SOURCE CODE

A.1. mainsim.h /* mainsim.h This header file includes general orthogonal cutting information and definitions. */ #define CHUNK_SIZE 256 #define STATES 4 #ifndef PI # define PI 3.1415928 #endif /* Defines the structure for tool information */ typedef struct { double mx, my, kx, ky, cx, cy, zetax, zetay; } TOOL_INFO; /* Defines the structure for workpiece information */ typedef struct { double N, d0, cutstiff; } WORKPIECE_INFO; /* Defines the structure for cutting parameters */ typedef struct { double DOC,b,theta,alpha,beta; } CUT_PARAM; /* Defines the structure to hold simulation parameters */ typedef struct { double hmax,hmin,h,epsmax,epsmin,tfi hmax,hmin,h,epsmax,epsmin,tfinal,V,feed,d0pi, nal,V,feed,d0pi, cx1,cx2,cx3,cy1,cy2,cy3; long int maxiter; } SIM_PARAM;

100 A.2. mainsim.c /* mainsim.c Thomas O'Reilly June 4 1997 Main file for orthogonal cutting simulation */ #include #include #include #include



#include "orthsim.h" #include "cutlist.h" /* Function prototypes */ void main(int argc, char *argv[]) { /* Declare variables */ double deltat,omegax,omegay; int revsteps, n, i, j; /* Declares new_history function */ extern HISTORY *first_history(); /* Declares main simulation function */ extern void orthsim(); /* Declares pointer to derivative function */ extern void orthder(); extern void stab_check1(); extern void stab_check2(); extern void stab_check3(); extern void stab_check4(); void (*p_schk)(); void (*p_ypfun)(); FILE *fp; HISTORY *p_hist; STEPDATA spres; clock_t elapsed; SIM_PARAM sim; TOOL_INFO tool; WORKPIECE_INFO work; CUT_PARAM cut; STAB_CHART_PARAM stab_p; double blim; /* Initialize variables */ elapsed=clock(); printf("main\n"); p_ypfun=orthder;

/* Function pointer for derivative function */

/* allocates space for and initializes pointer for first history structure */ p_hist=first_history(); init_step(&spres); /* The following initializes stuff */ /* Loads default parameters */

101 def_param(&sim, &tool, &work, &cut, &stab_p); /* Overrides defaults with contents of file simulate.cfg */ load_setup_file(&sim, &tool, &work, &cut, &stab_p); /* Overrides deafults or file values with command line arugments if they exist. */ if (argc>1) cmd_line_arg(argc,argv,&sim, cmd_line_arg(argc,argv,& sim, &tool, &work, &cut, &stab_p); /* Initializes parameters based upon values loaded at this point */ init_param(&sim, &tool, &work, &cut, &stab_p); /* Determines desired mode of program operation */ if (sim.mode==1) /* Calls the simulation function.*/ orthsim(p_ypfun,&sim,p_hist,&spres,0,tool,work,cut); else if (sim.mode==2) { /* Produce a stability chart for constant spindle speed, no tool motion */ p_schk=stab_check1; sweep_speeds(p_ypfun,p_schk,sim,&spres,p_hist,tool,work,cut,stab_p); } else if (sim.mode==3) /* Produce a stability chart based upon a maximum allowable level of vibration. Use for systems with forced tool motion */ { p_schk=stab_check2; sweep_speeds(p_ypfun,p_schk,sim,&spres,p_hist,tool,work,cut,stab_p); } else if (sim.mode==4) /* Produce a stability chart based upo Use for variable spindle speed and constant feed rate */ { p_schk=stab_check3; sweep_speeds(p_ypfun,p_schk,sim,&spres,p_hist,tool,work,cut,stab_p); } else if (sim.mode==5) /* Produce a stability chart with a minimum level of vibration Use for variable spindle speed and constant feed per revolution */ { p_schk=stab_check4; sweep_speeds(p_ypfun,p_schk,sim,&spres,p_hist,tool,work,cut,stab_p); } else if (sim.mode==6) /* Generates and saves the envelope of the chatter */ envsim(p_ypfun,sim,&spres,p_hist,tool,work,cut); /* Compute and display the elapsed time the program ran */ elapsed=clock(); printf("Elapsed Time: %f\n",(float) elapsed/CLOCKS_PER_SEC); /* Write the elapsed run time to a file */ /* fp=fopen("runtime.txt","w"); fprintf(fp,"Elapsed Time: %f\n",(float) elapsed/CLOCKS_PER_SEC); fclose(fp); */ exit(0); }

A.3. orthsim.h /* orthsim.h

102 This header file includes general orthogonal cutting information and definitions. */ /* Size of memory chunks */ #define CHUNK_SIZE 256 /* Number of state variables */ #define STATES 4 /* Number of state variable accelerations to track */ #define ACC_STATES 2 #ifndef PI # define PI 3.1415928 #endif /* Defines the structure for tool information */ typedef struct { double mx, my, kx, ky, cx, cy, zetax, zetay, omegax, omegay; } TOOL_INFO; /* Defines the structure for workpiece information */ typedef struct { double N, d0, cutstiff, rotperiod; } WORKPIECE_INFO; /* Defines the structure for cutting parameters */ typedef struct { double DOC,b,theta,alpha,beta,db_dt; } CUT_PARAM; /* Defines the structure to hold simulation parameters */ typedef struct { double hmax,hmin,h,epsmax,epsmin,tfi hmax,hmin,h,epsmax,epsmin,tfinal,V,feed,d0pi,int_time, nal,V,feed,d0pi,int_time, cx1,cx2,cx3,cy1,cy2,cy3,hinc,hdec,hmaxrat,hminrat,b,db_dt, rpmfreq,rpmamp,rpmcoeff,rpmomega,toolfreq,toolomega,toolamp, xmeandisp,xrange,xmindisp,xmaxdisp,feedcoeff,xmin,txmin,xmax,txmax, xvmin,txvmin,xvmax,txvmax,xamin,txamin,xamax,txamax,ymin,tymin,ymax, tymax,yvmin,tyvmin,yvmax,tyvmax,yamin,tyamin,yamax,tyamax; long int maxiter, totsteps; int int_cut, int_step, mode, fileskip, constfeed; } SIM_PARAM; /* Defines structure to hold stability chart parameters */ typedef struct { double Nmin,Nmax,bmin,bmax,bminabs,b Nmin,Nmax,bmin,bmax,bminabs,bmaxabs,increment,tolerance, maxabs,increment,tolerance, maxxmotion,slots; int Nsteps,maxchanges,initrot; } STAB_CHART_PARAM; /* Some #define #define #define #define #define #define #define

Number used to control printing and file handling */ PRINT 1 NO_PRINT 0 APPEND 1 NO_APPEND 0 WRITEFILE 1 NO_WRITEFILE 0 ROTATIONS 30

103 A.4. orthsim.c /* orthsim.c Thomas O'Reilly June 4 1997 Main file for orthogonal cutting simulation */ #include #include #include #include "orthsim.h" #include "cutlist.h" #include "integ.h" /* Function Prototypes */ void orthsim(void (*p_ypfun)(),SIM_PARAM *p_sim,HISTORY *p_hist, STEPDATA *p_spres,int nstart,TOOL_INFO tool,WORKPIECE_INFO work, CUT_PARAM cut); double orthpp(SIM_PARAM *p_sim, HISTORY *p_hpres,int n_pres,HISTORY **p_h1, int *n1, STEPDATA *p_spres,double Llast); void orthder(SIM_PARAM *p_sim, HISTORY *p_hpres,int n_pres,HISTORY **p_h1, int *n1,double Svecdot[STATES],STEPDATA *p_spres, double Llast); /* This function simulates Orthogonal Cutting */ void orthsim(void (*p_ypfun)(),SIM_PARAM *p_sim,HISTORY *p_hist, STEPDATA *p_spres,int nstart,TOOL_INFO tool,WORKPIECE_INFO work, CUT_PARAM cut) { /* Variable type declaration */ int nfinal, result; extern int rk45ad(); HISTORY p_hfinal; char *stepdatafile="stepdata.tx *stepdatafile="stepdata.txt"; t"; RK45CONST *p_rkcon;

/* Structure of rk45 constants */

/* Allocates space for rk45 constants */ p_rkcon=(RK45CONST *) malloc(sizeof(RK45CONST)); /* Initializes rk45 constants */ rk45init(p_rkcon); result=rk45ad(p_ypfun,p_rkcon,p_sim,p_hist,nstart,&p_hfinal, &nfinal,p_spres,PRINT,WRITEFILE,NO_APPEND,stepdatafile); if (result==0) { printf("\nThe simulation reached the maximum number of time steps\n"); printf("but didn't run to the desired final time.\n"); result=p_sim->maxiter; } else printf("\nThe simulation was successful.\n"); if (p_sim->int_cut!=0) { printf("The cut became interrupted at time step %d, Time = %f\n" ,p_sim->int_step,p_sim->int_time); }

104 return; } /* This function determines the thickness of the chip for orthogonal cutting simulations using linear interpolation. */ double orthpp(SIM_PARAM *p_sim, HISTORY *p_hpres,int npres,HISTORY **p_h1, int *p_n1, STEPDATA *p_spres,double Llast) /* The variables h_p1 and n_p1 refer to the first time step that bracketed the cutter position on the previous previous time step. These variables are updated in this function and the variables h_p2 and n_p2 are then set to refer to the previous position immediately on the other side of the tool during this time step */ { int n2,ind_temp,temp_ind,repeat n2,ind_temp,temp_ind,repeat; ; double n,Rdiff,toolmotion; HISTORY *p_h2,*hist_temp; toolmotion=0.0; n=0.0; if (p_sim->toolamp!=0.0) toolmotion=p_sim->toolamp*sin(2.0*PI*p_sim->toolfreq*p_spres->T); p_hpres->L[npres]=p_sim->V*p_spres->T+p_spres->Svec[2] +p_sim->rpmcoeff*sin(p_sim->rpmomega*p_spres->T) +toolmotion; p_hpres->R[npres]=p_sim->feed*p_spres->T+p_spres->Svec[0] +p_sim->feedcoeff*sin(p_sim->rpmomega*p_spres->T); /* This simulates a smooth entrance into the workpiece. */ if (p_hpres->L[npres]d0 (p_hpres->L[npres]d0pi) pi) { p_hpres->prev_cut[npres]=NULL; if (p_hpres->R[npres]>0.0) return p_hpres->R[npres]; else return 0.0; } else do { repeat=0; n=n+1.0; if ((**p_h1).L[*p_n1]>p_hpres->L ((**p_h1).L[*p_n1]>p_hpres->L[npres]-p_sim->d0pi*n) [npres]-p_sim->d0pi*n) { while ((**p_h1).L[*p_n1]>p_hpres-> ((**p_h1).L[*p_n1]>p_hpres->L[npres]-p_sim->d0pi*n) L[npres]-p_sim->d0pi*n) { *p_n1=*p_n1-1; if (*p_n1<0) { *p_n1=CHUNK_SIZE-1; *p_h1=(**p_h1).p_prev; } } } else { while ((**p_h1).L[*p_n1] ((**p_h1).L[*p_n1]L[npres]-p_sim->d0pi*n) L[npres]-p_sim->d0pi*n) { *p_n1=*p_n1+1; if (*p_n1>CHUNK_SIZE) { *p_n1=0;

105 *p_h1=(**p_h1).p_next; } } *p_n1=*p_n1-1; if (*p_n1<0) { *p_n1=CHUNK_SIZE-1; *p_h1=(**p_h1).p_prev; } } if (n==1) { hist_temp=*p_h1; ind_temp=*p_n1; } if ((**p_h1).prev_cut[*p_n1]!=NU ((**p_h1).prev_cut[*p_n1]!=NULL) LL) { temp_ind=*p_n1; *p_n1=(**p_h1).prev_ind[*p_n1]; *p_h1=(**p_h1).prev_cut[temp_ind]; repeat=1; } } while (repeat); n2=*p_n1+1; p_h2=*p_h1; if (n2>=CHUNK_SIZE) { n2=0; p_h2=(**p_h1).p_next; } Rdiff=p_hpres->R[npres]-((**p_h1).R[*p_n1]+ (p_hpres->L[npres]-p_sim->d0pi*n-(**p_h1).L[*p_n1])/ (p_h2->L[n2]-(**p_h1).L[*p_n1])*(p_h2->R[n2]-(**p_h1).R[*p_n1])); *p_h1=hist_temp; *p_n1=ind_temp; if (Rdiff<0.0) { Rdiff=0.0; p_hpres->prev_cut[npres]=hist_temp; p_hpres->prev_ind[npres]=ind_temp; } else { p_hpres->prev_cut[npres]=NULL; } return Rdiff; } /* This function determines the derivatives of the state variables for orthogonal cutting simulations. */ void orthder(SIM_PARAM *p_sim, HISTORY *p_hpres,int npres,HISTORY **p_h1, int *p_n1,double Svecdot[STATES],STEPDATA *p_spres, double Llast) { p_spres->RDIFF=orthpp(p_sim,p_hpres,npres,p_h1,p_n1,p_spres,Llast); /* Checks for interrupted cutting */ if (p_spres->RDIFF==0.0)

106 if ((p_hpres->R[npres]!=0.0)& ((p_hpres->R[npres]!=0.0)&&(p_sim->int_cut==0)) &(p_sim->int_cut==0)) { p_sim->int_cut=1; p_sim->int_time=p_spres->T; p_sim->int_step=p_sim->totsteps; } Svecdot[0]=p_spres->Svec[1]; Svecdot[1]=p_sim->cx1*p_spres->Svec[1] + p_sim->cx2*p_spres->Svec[0] + p_sim->cx3*(p_sim->b+p_simp_sim->cx3*(p_sim->b+p_sim->db_dt*p_spres->T)*p_spres->R >db_dt*p_spres->T)*p_spres->RDIFF; DIFF; Svecdot[2]=p_spres->Svec[3]; Svecdot[3]=p_sim->cy1*p_spres->Svec[3] + p_sim->cy2*p_spres->Svec[2] + p_sim->cy3*(p_sim->b+p_simp_sim->cy3*(p_sim->b+p_sim->db_dt*p_spres->T)*p_spres->R >db_dt*p_spres->T)*p_spres->RDIFF; DIFF; /* The next two lines keep track of accelerations */ p_spres->Svecdot[0]=Svecdot[1]; p_spres->Svecdot[1]=Svecdot[3]; return; }

A.5. cutlist.h /* cutlist.h This header file contains definitions of linked list structures*/ /* This structure defines the list elements used to record cutter path data*/ typedef struct history { double R[CHUNK_SIZE],L[CHUNK_SIZE]; struct history *p_next,*p_prev,*prev_cut *p_next,*p_prev,*prev_cut[CHUNK_SIZE]; [CHUNK_SIZE]; int prev_ind[CHUNK_SIZE]; } HISTORY; /* This defines the strucutre used to hold current position data. */ typedef struct stepdata { double Svec[STATES],Svecdot[STATES/2 Svec[STATES],Svecdot[STATES/2],T,RDIFF; ],T,RDIFF; } STEPDATA; /* This first structure is the definition for the surface list elements*/ typedef struct surface_patch { /* This defines the pointers to the previous and next list elements */ struct surface_patch *p_prev,*p_next; /* This defines the pointers to the list containing the start and end of the surface patch described by this element */ struct HISTORY *p_start,*p_end; /* This defines the indicies of the elements containing the position data */ int start_el,end_el; } SURFACE_PATCH;

107 A.6. cutlist.c /* cutlist.c This header file contains some of the functions for dealing with the linked lists. */ #include #include #include "orthsim.h" #include "cutlist.h" HISTORY *new_history(HISTORY *previous); HISTORY *first_history(void); void init_step(STEPDATA *p_spres); HISTORY *new_history(HISTORY *previous) /* This function creates a new turning history list element and attaches it to the previous list element and attaches the previous list element to the new element */ { HISTORY *p; int i,j; /* Allocates the space and exits if the space is not created. */ p=(HISTORY *) malloc(sizeof(HISTORY)); if (p==NULL) { printf("No new history history elements elements for you! Two megs.\n"); megs.\n"); exit(1); } /* Set the p_next pointer to point to nothing */ p->p_next=NULL; /* Set the p_prev pointer to point to the previous structure */ p->p_prev=previous; p->prev_cut[0]=NULL; /* The following block of code initializes all the elements in the history structure structure to zero. Since all the elements are filled before they are read this should be unnecessary */ /* for (i=0;iSvec[i][j]=0.0; p->Svecdot[i][j]=0.0; } p->L[i]=0.0; p->T[i]=0.0; p->R[i]=0.0; p->RDIFF[i]=0.0; } */ /* If this is not the first time step set the p_next pointer of the previous data structure to point to the newly produced structure */ if (previous!=NULL) previous->p_next=p; return p;

108 } HISTORY *first_history(void) /* This function creates the first history array and sets it up for all zero initial conditions */ { HISTORY *p; int j; p=new_history((HISTORY *) NULL); /* Allocates space for first history structure */ p->L[0]=0.0; p->R[0]=0.0; return p; } void init_step(STEPDATA *p_spres) /* Initializes the first time step when the program starts */ { int j; /* Set up initial conditions */ for (j=0;jSvec[j]=0.0; for (j=0;jSvecdot[j]=0.0; p_spres->T=0.0; p_spres->RDIFF=0.0; return; } SURFACE_PATCH *new_surface(SURFACE_PATCH *previous, SURFACE_PATCH *next) /* This function creates a new surface patch element and inserts it into the existing surface patch set. */ { SURFACE_PATCH *p; /* Allocates the space and exits if the space is not created. */ p=(SURFACE_PATCH *) malloc(sizeof(SURFACE_PAT malloc(sizeof(SURFACE_PATCH)); CH)); if (p==NULL) {printf("No new surface surface patches patches for you! Two megs.\n"); megs.\n"); exit(1);} p->p_next=next; p->p_prev=previous; if (previous!=NULL) previous->p_next=p; if (next!=NULL) next->p_prev=p; return p; }

A.7. input.h /* input.h This header file contains information related to reading simulation

109 parameters from a setup file and from the command line */ #define LINESIZE 100 #define NOFILE 0 #define READFILE 1

A.8. input.c /* input.c Thomas O'Reilly June 4 1997 functions to load parameters for orthogonal cutting simulation */ #include #include #include #include



#include "orthsim.h" #include "cutlist.h" #include "input.h" const char *SETUP_FILE="simulate.cfg"; /* Function prototypes */ void def_param(SIM_PARAM *p_sim, TOOL_INFO *p_tool, WORKPIECE_INFO *p_work, CUT_PARAM *p_cut, STAB_CHART_PARAM *p_stab_p); void set_param(char **pp_var, char **pp_val,SIM_PARAM *p_sim, TOOL_INFO *p_tool,WORKPIECE_INFO *p_work,CUT_PARAM *p_cut, STAB_CHART_PARAM *p_stab_p); void init_param(SIM_PARAM *p_sim, TOOL_INFO *p_tool, WORKPIECE_INFO *p_work, CUT_PARAM *p_cut, STAB_CHART_PARAM *p_stab_p); int load_setup_file(SIM_PARAM *p_sim, TOOL_INFO *p_tool, WORKPIECE_INFO *p_work, CUT_PARAM *p_cut, STAB_CHART_PARAM *p_stab_p); int cmd_line_arg(int argc, char *argv[], SIM_PARAM *p_sim, TOOL_INFO *p_tool, WORKPIECE_INFO *p_work, CUT_PARAM *p_cut, STAB_CHART_PARAM *p_stab_p); void make_token(char *string, char **pp_var, char **pp_val); void make_setup_file(SIM_PARAM *p_sim, TOOL_INFO *p_tool, WORKPIECE_INFO *p_work,CUT_PARAM *p_cut, STAB_CHART_PARAM *p_stab_p); /* Actual stuff starts here */ void def_param(SIM_PARAM *p_sim, TOOL_INFO *p_tool, WORKPIECE_INFO *p_work, CUT_PARAM *p_cut, STAB_CHART_PARAM *p_stab_p) /* Initializes variables using default parameters This function recieves pointers to the simulation, tool, workpiece, and cutting data structures to be initialized. */ { p_tool->kx=5500000.0; /* Tool stiffness in x direction */ p_tool->ky=3300000.0; /* Tool stiffness in y direction */ p_tool->mx=0.48; /* Effective mass in x direction */ p_tool->my=0.70; /* Effective mass in y direction */ p_tool->zetax=0.0216; /* Damping ratio for x direction */ p_tool->zetay=0.013; /* Damping ratio for y direction */

110 p_work->cutstiff=1e9; p_work->d0=0.0394; p_work->N=400.0;

/* Workpiece cutting stiffness /* Workpiece Diameter /* Rotational Speed, RPM

p_cut->DOC=0.0001; p_cut->b=0.001; p_cut->db_dt=0.0; p_cut->theta=PI/9.0;

/* /* /* /*

*/ */ */

Depth of cut, meters per revolution */ Width of cut */ Change in width of cut per unit time */ Angle of cutting force */

/* Defines some simulation parameters */ p_sim->int_cut=0; /* Checks for interrupted cutting */ p_sim->epsmin=0.0000001; /* lower bound on error before increasing step size*/ p_sim->epsmax=0.000001; /* Maximum error per step */ p_sim->maxiter=100000; /* Maximum number of time steps */ p_sim->totsteps=0; /* Counts total number of steps elapsed */ p_sim->tfinal=1.0; /* End time */ p_sim->hinc=1.5; /* Ratio by which step size is increased if error is too small */ p_sim->hdec=0.75; /* Ratio by which step size is decreaed if error is too large */ p_sim->hmaxrat=100.0; /* Ratio of maximum step size to initial step size */ p_sim->hminrat=100000.0; /* Ratio of mimimum step size to initial step size */ p_sim->constfeed=0; /* Determines whether the feed rate is constant or not*/ p_stab_p->bmax=0.01; /* Maximum initial width of cut */ p_stab_p->bmaxabs=0.02; /* Maximum final width of cut */ p_stab_p->bmin=0.0; /* Minimum initial width of cut */ p_stab_p->bminabs=0.0; /* Minimum final width of cut */ p_stab_p->slots=32.0; /* Increment to shift upper and lower limits by */ p_stab_p->tolerance=0.0000025; p_stab_p->tolerance=0.0000 025; /* Tolerance on width of cut */ p_stab_p->Nmax=10000.0; /* Maximum spindle speed */ p_stab_p->Nmin=1000.0; /* Minimum spindle speed */ p_stab_p->Nsteps=1000; /* Number of increments in spindle speed */ p_stab_p->initrot=5; /* Initial rotations */ p_stab_p->maxchanges=5; /* Number of changes in amplkitude to watch for */ p_stab_p->maxxmotion=1e-6; /* Maximum permissible x-direction motion amplitude*/ p_sim->rpmfreq=0.0; p_sim->rpmamp=0.0; p_sim->toolfreq=0.0; p_sim->toolamp=0.0;

/* /* /* /*

p_sim->fileskip=5;

/* Time step data will be saved every fileskip steps */ /* Mode of operation :1=time domain 2=stability type 1 3=stability type 2 4=envelope */

p_sim->mode=1;

Frequency of spindle speed variation */ Amplitude of spindle speed variation */ Frequency of forced tool y motion */ Aplitude of forced tool y motion */

return; } int load_setup_file(SIM_PARAM *p_sim, TOOL_INFO *p_tool, WORKPIECE_INFO *p_work, CUT_PARAM *p_cut,STAB_CHART_PARAM *p_stab_p) /* This function reads simulation parameters from a data file This function recieves pointers to the simulation, tool, workpiece, and cutting data structures to be initialized. */

111 { FILE *fp; /* string is the line read from the file, p_var is the varaible name string pulled out of that line, p_val is the value string pulled from that line. */ char string[LINESIZE],*p_var[LI string[LINESIZE],*p_var[LINESIZE],*p_val[LINESIZE]; NESIZE],*p_val[LINESIZE]; /* Opens input file for reading. */ fp=fopen(SETUP_FILE,"r"); /* Checks to make sure the file was read correctly. If not prints out a message and calls the function to create it. */ if (fp==NULL) { printf("Setup file %s not found. Creating it from defaults\n",SETUP_FILE); defaults\n",SETUP_FILE); printf("Simulation will run using default settings\n"); make_setup_file(p_sim, p_tool,p_work,p_cut,p_sta p_tool,p_work,p_cut,p_stab_p); b_p); return NOFILE; } /* Reads and scans through the input file. */ while (fgets(string,LINESIZE-1, (fgets(string,LINESIZE-1,fp)!=NULL) fp)!=NULL) { make_token(string,p_var,p_val); set_param(p_var,p_val,p_sim, set_param(p_var,p_val,p_ sim, p_tool,p_work,p_cut,p_stab_p p_tool,p_work,p_cut,p_stab_p); ); } fclose(fp); return READFILE; } int cmd_line_arg(int argc, char *argv[], SIM_PARAM *p_sim, TOOL_INFO *p_tool, WORKPIECE_INFO *p_work, CUT_PARAM *p_cut, STAB_CHART_PARAM *p_stab_p) /* This function function reads simulation parameters parameters from the the command line. argc, and argv are the command line variables passed to main. This function recieves pointers to the simulation, tool, workpiece, and cutting data structures to be initialized. */ { char *p_var[LINESIZE],*p_val[LI *p_var[LINESIZE],*p_val[LINESIZE]; NESIZE]; /* This loop analyzes each of the command line arguments */ while (--argc>0) { make_token(*++argv,p_var,p_val); set_param(p_var,p_val,p_sim, set_param(p_var,p_val,p_ sim, p_tool,p_work,p_cut,p_stab_p p_tool,p_work,p_cut,p_stab_p); ); } return 1; } void make_token(char *string, char **pp_var, char **pp_val) /* Splits the input variable 'string' into two parts, pp_var and pp_val using an '=' symbol to separate separate them. */ { /* This line pulls out the part of 'string' between the first non-space character and the first appearance of an '=' sign. */ *pp_var=strtok(string, " ="); /* This line pulls out the part of 'string' between the '=' sign and the next space. */ *pp_val=strtok(NULL," "); } void set_param(char **pp_var, char **pp_val,SIM_PARAM *p_sim,

112 TOOL_INFO *p_tool,WORKPIECE_INFO *p_work,CUT_PARAM *p_cut, STAB_CHART_PARAM *p_stab_p) /* Sets variables according to the variable name string referenced by pp_var, and the value string referenced by pp_val. This function recieves pointers to the simulation, tool, workpiece, and cutting data structures to be initialized. */ { if (!strcmp(*pp_var,"kx")) p_tool->kx=atof(*pp_val); if (!strcmp(*pp_var,"ky")) p_tool->ky=atof(*pp_val); if (!strcmp(*pp_var,"mx")) p_tool->mx=atof(*pp_val); if (!strcmp(*pp_var,"my")) p_tool->my=atof(*pp_val); if (!strcmp(*pp_var,"zetax")) p_tool->zetax=atof(*pp_val); p_tool->zetax=atof(*pp_val) ; if (!strcmp(*pp_var,"zetay")) p_tool->zetay=atof(*pp_val); p_tool->zetay=atof(*pp_val) ;

/* Tool stiffness in x direction */ /* Tool stiffness in y direction */ /* Effective mass in x direction */ /* Effective mass in y direction */ /* Damping ratio for x direction */ /* Damping ratio for y direction */

if (!strcmp(*pp_var,"cutstiff (!strcmp(*pp_var,"cutstiff")) ")) p_work->cutstiff=atof(*pp_val); p_work->cutstiff=atof(*pp_v al); /* Workpiece cutting stiffness */ if (!strcmp(*pp_var,"d0")) p_work->d0=atof(*pp_val); /* Workpiece Diameter */ if (!strcmp(*pp_var,"N")) p_work->N=atof(*pp_val); /* Rotational Speed, RPM */ if (!strcmp(*pp_var,"DOC")) p_cut->DOC=atof(*pp_val); if (!strcmp(*pp_var,"b")) p_cut->b=atof(*pp_val); if (!strcmp(*pp_var,"db_dt")) p_cut->db_dt=atof(*pp_val);

/* Depth of cut, meters per revolution */ /* Width of cut */ /* rate of change of Width of cut */

if (!strcmp(*pp_var,"theta")) p_cut->theta=atof(*pp_val)*PI/180.0; p_cut->theta=atof(*pp_val)* PI/180.0;

/* Angle of cutting cutting force */

/* Defines some simulation parameters */ if (!strcmp(*pp_var,"epsmin") (!strcmp(*pp_var,"epsmin")) ) p_sim->epsmin=atof(*pp_val); p_sim->epsmin=atof(*pp_val) ; /* lower bound on error before increasing step size*/ if (!strcmp(*pp_var,"epsmax") (!strcmp(*pp_var,"epsmax")) ) p_sim->epsmax=atof(*pp_val); p_sim->epsmax=atof(*pp_val) ; /* Maximum error per step */ if (!strcmp(*pp_var,"maxiter" (!strcmp(*pp_var,"maxiter")) )) p_sim->maxiter=atoi(*pp_val); p_sim->maxiter=atoi(*pp_val ); /* Maximum number of time steps */ if (!strcmp(*pp_var,"tfinal") (!strcmp(*pp_var,"tfinal")) ) p_sim->tfinal=atof(*pp_val); p_sim->tfinal=atof(*pp_val) ; /* End time */ if (!strcmp(*pp_var,"hinc")) p_sim->hinc=atof(*pp_val); /* Ratio by which step size is increased if error is too small */ if (!strcmp(*pp_var,"hdec")) p_sim->hdec=atof(*pp_val); /* Ratio by which step size is decreased if error is too large */ if (!strcmp(*pp_var,"hmaxrat" (!strcmp(*pp_var,"hmaxrat")) )) p_sim->hmaxrat=atof(*pp_val); p_sim->hmaxrat=atof(*pp_val ); /* Ratio of maximum step size to initial step size */ if (!strcmp(*pp_var,"hminrat" (!strcmp(*pp_var,"hminrat")) )) p_sim->hminrat=atof(*pp_val); p_sim->hminrat=atof(*pp_val ); /* Ratio of mimimum step size to initial step size */ if (!strcmp(*pp_var,"rpmfreq")) (!strcmp(*pp_var,"rpmfreq")) p_sim->rpmfreq=atof(*pp_val); if (!strcmp(*pp_var,"rpmamp") (!strcmp(*pp_var,"rpmamp")) )

/* Frequency Frequency of spindle speed variation */ /* Amplitude of spindle speed variation */

113 p_sim->rpmamp=atof(*pp_val); if (!strcmp(*pp_var,"toolfreq")) /* Frequency of forced tool y motion */ p_sim->toolfreq=atof(*pp_val); if (!strcmp(*pp_var,"toolamp")) (!strcmp(*pp_var,"toolamp")) /* Aplitude of forced tool y motion */ p_sim->toolamp=atof(*pp_val); if (!strcmp(*pp_var,"mode")) /* Mode of operation */ p_sim->mode=atoi(*pp_val); if (!strcmp(*pp_var,"fileskip (!strcmp(*pp_var,"fileskip")) ")) p_sim->fileskip=atoi(*pp_val); p_sim->fileskip=atoi(*pp_va l); /* step step saving saving interval */ if (!strcmp(*pp_var,"constfee (!strcmp(*pp_var,"constfeed")) d")) p_sim->constfeed=atoi(*pp_val); p_sim->constfeed=atoi(*pp_v al); /* Whether or not feed is constant */ if (!strcmp(*pp_var,"bmax")) p_stab_p->bmax=atof(*pp_val); p_stab_p->bmax=atof(*pp_val ); /* Maximum initial width of cut */ if (!strcmp(*pp_var,"bmaxabs" (!strcmp(*pp_var,"bmaxabs")) )) p_stab_p->bmaxabs=atof(*pp_val); p_stab_p->bmaxabs=atof(*pp_ val); /* Maximum final width of cut */ if (!strcmp(*pp_var,"bmin")) p_stab_p->bmin=atof(*pp_val); p_stab_p->bmin=atof(*pp_val ); /* Minimum initial width of cut */ if (!strcmp(*pp_var,"bminabs" (!strcmp(*pp_var,"bminabs")) )) p_stab_p->bminabs=atof(*pp_val); p_stab_p->bminabs=atof(*pp_ val); /* Minimum final width of cut */ if (!strcmp(*pp_var,"slots")) p_stab_p->slots=atof(*pp_val); p_stab_p->slots=atof(*pp_va l); /* Increment to shift upper and lower limits by */ if (!strcmp(*pp_var,"toleranc (!strcmp(*pp_var,"tolerance")) e")) p_stab_p->tolerance=atof(*pp_val); p_stab_p->tolerance=atof(*p p_val); /* Tolerance on width of cut */ if (!strcmp(*pp_var,"Nmax")) p_stab_p->Nmax=atof(*pp_val); p_stab_p->Nmax=atof(*pp_val ); /* Maximum spindle speed */ if (!strcmp(*pp_var,"Nmin")) p_stab_p->Nmin=atof(*pp_val); p_stab_p->Nmin=atof(*pp_val ); /* Minimum spindle speed */ if (!strcmp(*pp_var,"Nsteps") (!strcmp(*pp_var,"Nsteps")) ) p_stab_p->Nsteps=atoi(*pp_val); p_stab_p->Nsteps=atoi(*pp_v al); /* Number of increments in spindle speed */ if (!strcmp(*pp_var,"initrot" (!strcmp(*pp_var,"initrot")) )) p_stab_p->initrot=atoi(*pp_val); p_stab_p->initrot=atoi(*pp_ val); /* Initial rotations */ if (!strcmp(*pp_var,"maxchang (!strcmp(*pp_var,"maxchanges")) es")) p_stab_p->maxchanges=atoi(*pp_val); p_stab_p->maxchanges=atoi(* pp_val); /* Number of changes in amplitude to watch for */ if (!strcmp(*pp_var,"maxxmoti (!strcmp(*pp_var,"maxxmotion")) on")) p_stab_p->maxxmotion=atof(*pp_val); p_stab_p->maxxmotion=atof(* pp_val); /* Maximum permissible x-direction motion*/ return; } void init_param(SIM_PARAM *p_sim, TOOL_INFO *p_tool, WORKPIECE_INFO *p_work, CUT_PARAM *p_cut, STAB_CHART_PARAM *p_stab_p) /* Computes values based on the parameters defined. This function recieves pointers to the simulation, tool, workpiece, and cutting data structures to be initialized. */ { double deltat; /* Used as an intermediate variable for finding the time step */ /* Compute Damping Constant */ p_tool->cx=2.0*p_tool->zetax*sqrt(p_tool->mx*p_tool->kx); p_tool->cy=2.0*p_tool->zetay*sqrt(p_tool->my*p_tool->ky); p_work->rotperiod=60.0/p_work->N; p_work->rotperiod=60.0/p_w ork->N; p_cut->alpha=sin(p_cut->theta); p_cut->alpha=sin(p_cut->theta); p_cut->beta=cos(p_cut->theta); p_cut->beta=cos(p_cut->the ta);

/* Rotation period */

/* X-direction force component /* Y-direction force component

*/ */

/* Compute constants used in derivative function. If the stiffness or mass in any direction is zero it is assumed that that DOF is inactive. */

114 if ((p_tool->kx==0.0)||(p_tool((p_tool->kx==0.0)||(p_tool->mx==0.0)) >mx==0.0)) { p_sim->cx1=0.0; p_sim->cx2=0.0; p_sim->cx3=0.0; p_tool->omegax=0.0; } else { p_sim->cx1=-p_tool->cx/p_tool->mx; p_sim->cx2=-p_tool->kx/p_tool->mx; /* This is not actually cx3 but is multiplied by p_cut->b to get cx3 */ p_sim->cx3=-p_cut->alpha*p_work->cutstiff/p_tool->mx; p_tool->omegax=sqrt(p_tool->kx/p_tool->mx); } if ((p_tool->ky==0.0)||(p_tool((p_tool->ky==0.0)||(p_tool->my==0.0)) >my==0.0)) { p_sim->cy1=0.0; p_sim->cy2=0.0; p_sim->cy3=0.0; p_tool->omegay=0.0; } else { p_sim->cy1=-p_tool->cy/p_tool->my; p_sim->cy2=-p_tool->ky/p_tool->my; /* This is not actually cy3 but is multiplied by p_cut->b to get cy3 */ p_sim->cy3=-p_cut->beta*p_work->cutstiff/p_tool->my; p_tool->omegay=sqrt(p_tool->my/p_tool->my); } /* The next few lines of code sets the initial time step to be 1/10 the period of the highest natural freqency */ if ((p_tool->omegax==0.0)&&(p_t ((p_tool->omegax==0.0)&&(p_tool->omegay==0.0)) ool->omegay==0.0)) deltat=p_sim->tfinal/((double) deltat=p_sim->tfinal/((d ouble) p_sim->maxiter); else { if (p_tool->omegax>p_tool->om (p_tool->omegax>p_tool->omegay) egay) deltat=2.0*PI/p_tool->omegax/10.0; else deltat=2.0*PI/p_tool->omegay/10.0; } /* Defines some simulation parameters */ p_sim->h=deltat; p_sim->hmax=deltat*p_sim->hmaxrat; p_sim->hmin=deltat/p_sim->hminrat; /* Copies some things into p_sim so that fewer variables are passed around later */ /* Width of cut stuff */ p_sim->b=p_cut->b; p_sim->db_dt=p_cut->db_dt; /* workpiece diameter * pi. Used to find find last DOC. */ p_sim->d0pi=p_work->d0*PI; /* Linear Feed speed into workpiece */ p_sim->feed=p_cut->DOC*p_work->N/60.0; /* Linear relative velocity of tool wrt workpiece */ p_sim->V=p_work->N/60.0*p_work->d0*PI; /* Forced tool motion frequency in rad/sec */

115 p_sim->toolomega=p_sim->toolfreq*2.0*PI; /* Variation of spindle speed coefficient for computing L */ if (p_sim->rpmfreq!=0.0) { p_sim->rpmomega=p_sim->rpmfreq*2.0*PI; p_sim->rpmcoeff=p_sim->rpmamp*p_sim->V/p_sim->rpmomega; if (p_sim->constfeed==1) p_sim->feedcoeff=0.0; else p_sim->feedcoeff=p_sim->rpmamp*p_sim->feed/p_sim->rpmomega; } else { p_sim->rpmomega=0.0; p_sim->rpmcoeff=0.0; p_sim->feedcoeff=0.0; } /* Spindle speed increment for stability chart generation */ p_stab_p->increment=(p_stab_p->bmaxabs-p_stab_p->bmina p_stab_p->increment=(p_sta b_p->bmaxabs-p_stab_p->bminabs)/p_stab_p->slots; bs)/p_stab_p->slots; if (p_tool->kx!=0.0) p_sim->xmeandisp=-p_sim->b*p_work->cutstiff*p_cut->DOC/p_tool->kx *sin(p_cut->theta); else p_sim->xmeandisp=0.0; p_sim->xmaxdisp=p_sim->xmeandisp/(1+p_sim->rpmamp); p_sim->xmindisp=p_sim->xmeandisp/(1-p_sim->rpmamp); return; } void make_setup_file(SIM_PARAM *p_sim, TOOL_INFO *p_tool, WORKPIECE_INFO *p_work,CUT_PARAM *p_cut, STAB_CHART_PARAM *p_stab_p) /* Creates a configuration file if one is not found. Initializes variables using default parameters */ { FILE *fp; fp=fopen(SETUP_FILE,"w"); fprintf(fp,"kx=%e /* Tool stiffness in x direction */\n",p_tool->kx); fprintf(fp,"ky=%e /* Tool stiffness in y direction */\n",p_tool->ky); fprintf(fp,"mx=%f /* Effective mass in x direction */\n",p_tool->mx); fprintf(fp,"my=%f /* Effective mass in y direction */\n",p_tool->my); fprintf(fp,"zetax=%e /* Damping ratio for x direction */\n", p_tool->zetax); fprintf(fp,"zetay=%e /* Damping ratio for y direction */\n", p_tool->zetay); fprintf(fp,"cutstiff=%e /* Workpiece cutting stiffness */\n", p_work->cutstiff); fprintf(fp,"d0=%f /* Workpiece Diameter */\n",p_work->d0); fprintf(fp,"N=%f /* Rotational Speed, RPM */\n",p_work->N); fprintf(fp,"DOC=%e /* Depth of cut, meters per revolution */\n", p_cut->DOC); fprintf(fp,"b=%e /* Width of cut */\n",p_cut->b); fprintf(fp,"db_dt=%e /* Rate of of change of width of of cut */\n",p_cut->db_dt); */\n",p_cut->db_dt); fprintf(fp,"theta=%f /* Angle of cutting force in degrees* /\n", p_cut->theta*180/PI);

116 /* Defines some simulation parameters */ fprintf(fp,"epsmin=%e /* lower bound on error before increasing\n", p_sim->epsmin); fprintf(fp,"epsmax=%e /* Maximum error per step */\n",p_sim->epsmax); fprintf(fp,"maxiter=%d /* Maximum number of time steps */\n", p_sim->maxiter); fprintf(fp,"tfinal=%f /* End time */\n",p_sim->tfinal); fprintf(fp,"hinc=%f /* Ratio by which step size is increased if\n", p_sim->hinc); fprintf(fp,"hdec=%f /* Ratio by which step size is decreaed if\n", p_sim->hdec); fprintf(fp,"hmaxrat=%f /* Ratio of maximum step size to initial\n", p_sim->hmaxrat); fprintf(fp,"hminrat=%f /* Ratio of mimimum step size to initial\n", p_sim->hminrat); fprintf(fp,"rpmfreq=%f /* Frequency of spindle speed variation variation */\n", p_sim->rpmfreq); fprintf(fp,"rpmamp=%f /* Amplitude of spindle speed variation */\n", p_sim->rpmamp); fprintf(fp,"toolfreq=%f /* Frequency of forced tool y motion */\n", p_sim->toolfreq); fprintf(fp,"toolamp=%f /* Aplitude of forced tool y motion motion */ \n", \n", p_sim->toolamp); fprintf(fp,"mode=%f /* Mode Mode of operation */ \n", p_sim->mode); fprintf(fp,"fileskip=%f /* Step Step save interval */ \n", p_sim->fileskip); fprintf(fp,"constfeed=%f /* Whther feed rate is constant \n", p_sim->constfeed); fprintf(fp,"bmax=%f /* Maximum initial width of cut \n", p_stab_p->bmax); fprintf(fp,"bmaxabs=%f /* Maximum final width of cut \n", p_stab_p->bmaxabs); fprintf(fp,"bmin=%f; /* Minimum initial width of cut \n", p_stab_p->bmin); fprintf(fp,"bminabs=%f /* Minimum final width of cut \n", p_stab_p->bminabs); fprintf(fp,"slots=%f /* Increment to shift upper and lower limits by \n", p_stab_p->slots); fprintf(fp,"tolerance=%f /* Tolerance on width of cut \n", p_stab_p->tolerance); fprintf(fp,"Nmax=%f /* Maximum spindle speed \n", p_stab_p->Nmax); fprintf(fp,"Nmin=%f /* Minimum spindle speed \n", p_stab_p->Nmin); fprintf(fp,"Nsteps=%f /* Number of increments in spindle speed \n", p_stab_p->Nsteps); fprintf(fp,"initrot=%f /* Initial rotations */\n", p_stab_p->initrot); fprintf(fp,"maxchanges=%f /* Number of changes in amplitude to watch for */\n", p_stab_p->maxchanges); fprintf(fp,"maxxmotion=%f /* Maximum permissible x-direction motion*/\n", p_stab_p->maxxmotion); fclose(fp); return; }

117 A.9. integ.h /* integ.h This header file contains information for the RK45 integrator function */ /* The following structure is used to hold constants for the rk45 function */ typedef struct { double c20,c21,c30,c31,c32,c40,c41,c c20,c21,c30,c31,c32,c40,c41,c42,c43,c51,c52,c53,c54, 42,c43,c51,c52,c53,c54, c60,c61,c62,c63,c64,c65,a1,a2,a3,a4,a5,b1,b2,b3,b4,b5,b6; } RK45CONST;

A.10. integ.c /* integ.c This file contains RK45 related functions for time domain simulation of IVP differential equations. */ #include #include #include #include "orthsim.h" #include "cutlist.h" #include "integ.h" /* Function prototypes */ void rk45init(RK45CONST *p_rkcon); double rk45step(void *p_ypfun(),RK45CONST *p_rkcon,SIM_PARAM *p_sim, int npres,int nnext,HISTORY *p_hpres,HISTORY *p_hnext, HISTORY **pp_hlastcut,int *p_nlastcut, STEPDATA *p_spres,STEPDATA *p_snext); int rk45ad(void *p_ypfun(),RK45CONST *p_rkcon,SIM_PARAM *p_sim, HISTORY *p_hpres,int npres,HISTORY **pp_hfinal,int *p_nfinal, STEPDATA *p_spres, int printout,int writefile,int append, char *datafile); void rk45init(RK45CONST *p_rkcon) /* This function initializes the RK45 coefficients */ { p_rkcon->c20=0.25; p_rkcon->c21=0.25; p_rkcon->c30=0.375; p_rkcon->c31=0.09375; p_rkcon->c32=0.28125; p_rkcon->c40=12.0/13.0; p_rkcon->c41=1932.0/2197.0; p_rkcon->c42=-7200.0/2197.0; p_rkcon->c43=7296.0/2197.0; p_rkcon->c51=439.0/216.0; p_rkcon->c52=-8.0; p_rkcon->c53=3680.0/513.0;

118 p_rkcon->c54=-845.0/4104.0; p_rkcon->c60=0.5; p_rkcon->c61=-8.0/27.0; p_rkcon->c62=2.0; p_rkcon->c63=-3544.0/2565.0; p_rkcon->c64=1859.0/4104.0; p_rkcon->c65=-0.275; p_rkcon->a1=25.0/216.0; p_rkcon->a2=0.0; p_rkcon->a3=1408.0/2565.0; p_rkcon->a4=2197.0/4104.0; p_rkcon->a5=-0.2; p_rkcon->b1=16.0/135.0; p_rkcon->b2=0.0; p_rkcon->b3=6656.0/12825.0; p_rkcon->b4=28561.0/56430.0; p_rkcon->b5=-0.18; p_rkcon->b6=2.0/55.0; return; } double rk45step(void *p_ypfun(),RK45CONST *p_rkcon,SIM_PARAM *p_sim, int npres,int nnext,HISTORY *p_hpres,HISTORY *p_hnext, HISTORY **pp_hlastcut,int *p_nlastcut, STEPDATA *p_spres,STEPDATA *p_snext) /* This function makes a single step using RK45 and returns an estimate of the error involved in taking that step. It recieves as arguments: p_ypfun: pointer to derivative function p_rkcon: constants for rk45 functions p_sim: pointer to structure with simulation parameters npres: the index within a history structure of the current time step nnext: the index within a history structure of the next time step p_hpres: pointer to current history structure p_hnext: pointer to next history structure pp_hlastcut: pointer to the history structure containing the surface data that was involved in chip calculations the previous time step p_nlastcut: the most recent time step within pp_hlastcut that that is older older than the true previous position on the previous time step p_spres: Present time step data */ { double F1[STATES],F2[STATES],F3[STAT F1[STATES],F2[STATES],F3[STATES],F4[STATES],F5[STATES],F6 ES],F4[STATES],F5[STATES],F6[STATES], [STATES], h4[STATES]; /* various various intermediate intermediate derivative derivative terms and 4th order estimate of next position */ double eps; /* error estimate */ double htemp; /* Temporary variable to hold time step size */ double Llast; /* Value of L at the start of the time step */ int i; /* used for counting in loops */ eps=0.0; /* Start of RK45 method code */ F1[0]=p_spres->Svec[1]*p_sim->h; F1[1]=p_spres->Svecdot[0]*p_sim->h; F1[2]=p_spres->Svec[3]*p_sim->h; F1[3]=p_spres->Svecdot[1]*p_sim->h; htemp=p_sim->h; Llast=p_hpres->L[npres];

119 p_sim->h=htemp*p_rkcon->c20; p_snext->T=p_spres->T+p_sim->h; for (i=0;iSvec[i]=p_spres->Svec[i]+p_rkcon->c21*F1[i]; p_ypfun(p_sim,p_hnext,nnext,pp_hlastcut,p_nlastcut,F2,p_snext,Llast); for (i=0;ih=htemp*p_rkcon->c30; p_snext->T=p_spres->T+p_sim->h; for (i=0;iSvec[i]=p_spres->Svec[i]+p_rkcon->c31*F1[i] +p_rkcon->c32*F2[i]; p_ypfun(p_sim,p_hnext,nnext,pp_hlastcut,p_nlastcut,F3,p_snext,Llast); for (i=0;ih=htemp*p_rkcon->c40; p_snext->T=p_spres->T+p_sim->h; for (i=0;iSvec[i]=p_spres->Svec[i]+p_rkcon->c41*F1[i] +p_rkcon->c42*F2[i]+p_rkcon->c43*F3[i]; p_ypfun(p_sim,p_hnext,nnext,pp_hlastcut,p_nlastcut,F4,p_snext,Llast); for (i=0;ih=htemp; p_snext->T=p_spres->T+p_sim->h; for (i=0;iSvec[i]=p_spres->Svec[i]+p_rkcon->c51*F1[i] +p_rkcon->c52*F2[i]+p_rkcon->c53*F3[i] +p_rkcon->c54*F4[i]; p_ypfun(p_sim,p_hnext,nnext,pp_hlastcut,p_nlastcut,F5,p_snext,Llast); for (i=0;ih=htemp*p_rkcon->c60; p_snext->T=p_spres->T+p_sim->h; for (i=0;iSvec[i]=p_spres->Svec[i]+p_rkcon->c61*F1[i] +p_rkcon->c62*F2[i]+p_rkcon->c63*F3[i] +p_rkcon->c64*F4[i]+p_rkcon->c65*F5[i]; p_ypfun(p_sim,p_hnext,nnext,pp_hlastcut,p_nlastcut,F6,p_snext,Llast); for (i=0;ih=htemp; for (i=0;iSvec[i]+p_rkcon->a1*F1[i]+p_rkcon->a3*F3[i] +p_rkcon->a4*F4[i]+p_rkcon->a5*F5[i]; p_snext->Svec[i]=p_spres->Svec[i]+p_rkcon->b1*F1[i] +p_rkcon->b3*F3[i]+p_rkcon->b4*F4[i] +p_rkcon->b5*F5[i]+p_rkcon->b6*F6[i]; eps=eps+fabs(h4[i]-p_snext->Svec[i]); } p_snext->T=p_spres->T+p_sim->h; p_ypfun(p_sim,p_hnext,nnext,pp_hlastcut,p_nlastcut,F1,p_snext,Llast); return eps; }

120 int rk45ad(void *p_ypfun(),RK45CONST *p_rkcon,SIM_PARAM *p_sim, HISTORY *p_hpres,int npres,HISTORY **pp_hfinal,int *p_nfinal, STEPDATA *p_spres, int printout,int writefile,int append, char *datafile) /* This is the function that uses an adaptive RK45 scheme to solve an IVP */ /* p_ypfun: pointer to derivative function p_sim: pointer to structure containing simulation parameters p_hpres: pointer to initial data structure for holding history info npres: index of initial time step within p_hpres pp_hfinal: points to the structure containing the last time step performed p_nfinal: index of last step performed within pp_hfinal pp_hfinal printout: if != 0 will print the progress of the run during execution writefile: if !=0 will open and write data to a file append: if !=0 will append to datafile datafile: contains name of output output file */ { int nnext,nlastcut,j; /* nnext: next time step index nlastcut: index of tool previous position on previous time step */ double htemp; /* Used to store the last regular step before the final step */ HISTORY *p_hnext,*p_hlastcut;

STEPDATA snext,*p_snext; STEPDATA *p_stemp,*p_sorig; double eps;

/* p_hnext: pointer to history structure for next time step p_hlastcut: pointer to structure for nlastcut. */ /* pointers to data for next time step */ /* Error estimate */

extern HISTORY *new_history(HISTORY *); /* new history structure function */ FILE *fp; if (writefile) { if (append==0) { fp=fopen(datafile, "w"); fprintf(fp,"%e, ",p_spres->T); for (j=0;jSvec[j]); fprintf(fp,"%10.9e, ", p_spres->Svecdot[0]); fprintf(fp,"%10.9e, ", p_spres->Svecdot[1]); fprintf(fp,"%10.9e, ", p_spres->RDIFF); fprintf(fp,"%10.9e, ",p_hpres->L[npres]); fprintf(fp,"%10.9e\n",p_hpres->R[npres]); } else fp=fopen(datafile, "a"); if (fp==NULL) printf("Error opening output file %s\n",datafile); } /* Keeps track of the original step data structure */ p_sorig=p_spres; p_snext=&snext; eps=0.0; /* Initializes the variables for the next time step */

121 p_hnext=p_hpres; nnext=npres+1; /* Produces a new history structure as needed */ if (nnext==CHUNK_SIZE) { if (p_hpres->p_next==NULL) p_hnext=(HISTORY *) new_history(p_hpres); else p_hnext=p_hpres->p_next; nnext=0; } /* Initializes variables for last cut data */ p_hlastcut=p_hpres; nlastcut=npres; /* Open data file for writing. the same name */

Will overwrite any existing file with

/* Loops until the final time or maximum iterations is reached */ while ((p_spres->Ttfina ((p_spres->Ttfinal)&&(p_sim->totsteps l)&&(p_sim->totstepsmaxiter)) maxiter)) { /* Corrects step sizes */ if (p_sim->h>p_sim->hmax) p_sim->h=p_sim->hmax; if (p_sim->hhmin) p_sim->h=p_sim->hmin; if ((p_sim->h+p_spres->T)>p_s ((p_sim->h+p_spres->T)>p_sim->tfinal) im->tfinal) { htemp=p_sim->h; p_sim->h=p_sim->tfinal-p_spres->T; } /* Produces a new history structure as needed */ if (nnext==CHUNK_SIZE) { if (p_hpres->p_next==NULL) { p_hnext=(HISTORY *) new_history(p_hpres); } else p_hnext=p_hpres->p_next; nnext=0; } /* printf("Time = %e, Step = %d, Step size = %e, ",p_hpres>T[npres],totsteps,p_sim->h);*/ eps=rk45step(p_ypfun,p_rkcon,p_sim,npres,nnext,p_hpres,p_hnext, &p_hlastcut,&nlastcut,p_spres,p_snext); if (printout) printf("%4.1f%%, Totstep= %d\r", p_spres->T/p_sim->tfinal*10 p_spres->T/p_sim->tfinal*100, 0, p_sim->totsteps); /* Updates everything to the next time step if the error is acceptable*/ if ((eps<=p_sim->epsmax)||(p_ ((eps<=p_sim->epsmax)||(p_sim->h<=p_sim->hmin)) sim->h<=p_sim->hmin)) { p_sim->totsteps=p_sim->totsteps+1; /* Points to the next present history data step */ npres=nnext; p_hpres=p_hnext; nnext=nnext+1;

122 /* Swaps pointers to the step data */ p_stemp=p_spres; p_spres=p_snext; p_snext=p_stemp; /* This write the output data when needed */ if (!(p_sim->totsteps%p_sim->fil (!(p_sim->totsteps%p_sim->fileskip)&&(writefile)) eskip)&&(writefile)) { fprintf(fp,"%e, ",p_spres->T); for (j=0;jSvec[j]); fprintf(fp,"%10.9e, ", p_spres->Svecdot[0]); fprintf(fp,"%10.9e, ", p_spres->Svecdot[1]); fprintf(fp,"%10.9e, ", p_spres->RDIFF); fprintf(fp,"%10.9e, ",p_hpres->L[npres]); fprintf(fp,"%10.9e\n",p_hpres->R[npres]); } if (p_spres->Svec[0]>(p_sim->xma (p_spres->Svec[0]>(p_sim->xmax)) x)) { p_sim->xmax=p_spres->Svec[0]; p_sim->txmax=p_spres->T; } if (p_spres->Svec[0]<(p_sim->xmi (p_spres->Svec[0]<(p_sim->xmin)) n)) { p_sim->xmin=p_spres->Svec[0]; p_sim->txmin=p_spres->T; } if (p_spres->Svec[1]>(p_sim->xvm (p_spres->Svec[1]>(p_sim->xvmax)) ax)) { p_sim->xvmax=p_spres->Svec[1]; p_sim->txvmax=p_spres->T; } if (p_spres->Svec[1]<(p_sim->xvm (p_spres->Svec[1]<(p_sim->xvmin)) in)) { p_sim->xvmin=p_spres->Svec[1]; p_sim->txvmin=p_spres->T; } if (p_spres->Svecdot[0]>(p_sim-> (p_spres->Svecdot[0]>(p_sim->xamax)) xamax)) { p_sim->xamax=p_spres->Svecdot[0]; p_sim->txamax=p_spres->T; } if (p_spres->Svecdot[0]<(p_sim-> (p_spres->Svecdot[0]<(p_sim->xamin)) xamin)) { p_sim->xamin=p_spres->Svecdot[0]; p_sim->txamin=p_spres->T; } if (p_spres->Svec[2]>(p_sim->yma (p_spres->Svec[2]>(p_sim->ymax)) x)) { p_sim->ymax=p_spres->Svec[2]; p_sim->tymax=p_spres->T; } if (p_spres->Svec[2]<(p_sim->ymi (p_spres->Svec[2]<(p_sim->ymin)) n)) { p_sim->ymin=p_spres->Svec[2]; p_sim->tymin=p_spres->T; } if (p_spres->Svec[3]>(p_sim->yvm (p_spres->Svec[3]>(p_sim->yvmax)) ax)) { p_sim->yvmax=p_spres->Svec[3]; p_sim->tyvmax=p_spres->T; } if (p_spres->Svec[3]<(p_sim->yvm (p_spres->Svec[3]<(p_sim->yvmin)) in))

123 { p_sim->yvmin=p_spres->Svec[3]; p_sim->tyvmin=p_spres->T; } if (p_spres->Svecdot[1]>(p_sim-> (p_spres->Svecdot[1]>(p_sim->yamax)) yamax)) { p_sim->yamax=p_spres->Svecdot[1]; p_sim->tyamax=p_spres->T; } if (p_spres->Svecdot[1]<(p_sim-> (p_spres->Svecdot[1]<(p_sim->yamin)) yamin)) { p_sim->yamin=p_spres->Svecdot[1]; p_sim->tyamin=p_spres->T; } } /* Increases the step size if the error estimate is too small */ if (epsepsmin) p_sim->h=p_sim->h*p_sim->hinc; /* Decreases the step size if the step is not acceptable */ if (eps>p_sim->epsmax) { p_sim->h=p_sim->h*p_sim->hdec; } } /* Sets the pointers to the final time step performed */ *pp_hfinal=p_hpres; *p_nfinal=npres; /* resets time step to the stored temporary value */ p_sim->h=htemp; /* Resets the original present time step values to the actual last time step*/ if (p_spres!=p_sorig) { p_sorig->T=p_spres->T; p_sorig->RDIFF=p_spres->RDIFF; p_sorig->Svecdot[0]=p_spres->Svecdot[0]; p_sorig->Svecdot[1]=p_spres->Svecdot[1]; p_sorig->Svec[0]=p_spres->Svec[0]; p_sorig->Svec[1]=p_spres->Svec[1]; p_sorig->Svec[2]=p_spres->Svec[2]; p_sorig->Svec[3]=p_spres->Svec[3]; } /* Close Data file */ if (writefile) fclose(fp); /* Checks if the simulation ran to the desired time */ if (p_spres->Ttfinal) return 0; else return p_sim->totsteps; }

124 A.11. stabchk.h /* stabchk.h Thomas O'Reilly June 4 1997 Stability checking stuff */ #define STABLE -1 #define UNSTABLE 1 #define STAB_ERROR 0

/* Value returned for stable system */ /* Value returned for unstable system */ /* Value returned if error detected */

A.12. stabchk.c /* stabchk.c Thomas O'Reilly June 4 1997 Stability checking stuff */ #include #include #include #include #include #include #include

"orthsim.h" "cutlist.h" "stabchk.h" "integ.h"

/* Function Prototypes */ int stab_check1(void (*p_ypfun)(),RK45CONST *p_rkcon,SIM_PARAM sim, HISTORY *p_hist,STEPDATA *p_spres, WORKPIECE_INFO *p_work, int printout,STAB_CHART_PARAM *p_stabp); int stab_check2(void (*p_ypfun)(),RK45CONST *p_rkcon,SIM_PARAM sim, HISTORY *p_hist,STEPDATA *p_spres, WORKPIECE_INFO *p_work, int printout,STAB_CHART_PARAM *p_stabp); int stab_check3(void (*p_ypfun)(),RK45CONST *p_rkcon,SIM_PARAM sim, HISTORY *p_hist,STEPDATA *p_spres, WORKPIECE_INFO *p_work, int printout,STAB_CHART_PARAM *p_stabp); int stab_check4(void (*p_ypfun)(),RK45CONST *p_rkcon,SIM_PARAM sim, HISTORY *p_hist,STEPDATA *p_spres, WORKPIECE_INFO *p_work, int printout,STAB_CHART_PARAM *p_stabp); double stab_limit(void (*p_ypfun)(),int (*p_schk)(),RK45CONST *p_rkcon, SIM_PARAM sim, HISTORY *p_hist,TOOL_INFO tool,WORKPIECE_INFO work, CUT_PARAM cut,STAB_CHART_PARAM *stab_p, STEPDATA *p_spres); void sweep_speeds(void (*p_ypfun)(),int (*p_schk)(),SIM_PARAM sim, STEPDATA *p_spres, HISTORY *p_hist,TOOL_INFO tool,WORKPIECE_INFO work, CUT_PARAM cut,STAB_CHART_PARAM stab_p); /* Real stuff starts here */ void sweep_speeds(void (*p_ypfun)(),int (*p_schk)(),SIM_PARAM sim, STEPDATA *p_spres, HISTORY *p_hist,TOOL_INFO tool,WORKPIECE_INFO work, CUT_PARAM cut,STAB_CHART_PARAM stab_p)

125 /* This function scans across a number of spindle speeds and examines the stability of the system at each. */ { double Ninc,blim; int Nsteps; FILE *fp; int i; RK45CONST *p_rkcon;

/* Structure of rk45 constants */

/* Allocates space for rk45 constants */ p_rkcon=(RK45CONST *) malloc(sizeof(RK45CONST)); /* Initializes rk45 constants */ rk45init(p_rkcon); /* Open Output file */ fp=fopen("stabch.txt","w"); if (fp==NULL) printf("stabch.txt file not opened"); work.N=stab_p.Nmin; Ninc=(double) (stab_p.Nmax-stab_p.Nmin)/ (stab_p.Nmax-stab_p.Nmin)/stab_p.Nsteps; stab_p.Nsteps; init_param(&sim, &tool, &work, &cut,&stab_p); printf("N=%f\n",work.N); while (work.N<=stab_p.Nmax) { blim=stab_limit(p_ypfun,p_schk,p_rkcon,sim,p_hist,to blim=stab_limit(p_ypfun, p_schk,p_rkcon,sim,p_hist,tool,work,cut,&stab_p, ol,work,cut,&stab_p, p_spres); fprintf(fp,"%e %e\n",work.N,blim); printf("%4.1f %% done, N=%f, blim=%f\n", (work.N-stab_p.Nmin)/(stab_p.Nmax-stab_p.Nmin)*100.0,work.N,blim); work.N=work.N+Ninc; init_param(&sim, &tool, &work, &cut,&stab_p); } fclose(fp); } double stab_limit(void (*p_ypfun)(),int (*p_schk)(),RK45CONST *p_rkcon,SIM_PARAM sim, HISTORY *p_hist,TOOL_INFO tool,WORKPIECE_INFO work, CUT_PARAM cut,STAB_CHART_PARAM *p_stabp, STEPDATA *p_spres) /* This Function finds the maximum value of b for stable cutting for the given system at a certain speed */ { int maxunstab,minunstab,meanuns maxunstab,minunstab,meanunstab; tab; double range,bmean; cut.b=p_stabp->bmax; sim.b=p_stabp->bmax; init_param(&sim, &tool, &work, &cut,p_stabp); maxunstab=p_schk(p_ypfun,p_rkcon,sim,p_hist,p_spres,&work, NO_PRINT,p_stabp); cut.b=p_stabp->bmin; sim.b=p_stabp->bmin; init_param(&sim, &tool, &work, &cut,p_stabp); minunstab=p_schk(p_ypfun,p_rkcon,sim,p_hist,p_spres,&work, NO_PRINT,p_stabp);

126 while (maxunstab==minunstab) { /* Check to make sure a limit is bracketed */ if ((p_stabp->bmax==p_stabp-> ((p_stabp->bmax==p_stabp->bmaxabs)&&(p_stabp->bmin==p_ bmaxabs)&&(p_stabp->bmin==p_stabp->bminabs)) stabp->bminabs)) { printf("No limit bracketed. Thats my story and I'm sticking by it.\n"); return p_stabp->bmax; } /* If both upper and lower limits are stable, raise both */ else if (maxunstab==STABLE) { p_stabp->bmin=p_stabp->bmax; p_stabp->bmax=p_stabp->bmax+p_stabp->increment; cut.b=p_stabp->bmax; sim.b=p_stabp->bmax; if (cut.b>p_stabp->bmaxabs) { cut.b=p_stabp->bmaxabs; sim.b=p_stabp->bmaxabs; } init_param(&sim, &tool, &work, &cut,p_stabp); maxunstab=p_schk(p_ypfun,p_rkcon,sim,p_hist,p_spres,&work, NO_PRINT,p_stabp); } /* If both upper and lower limits are unstable, lower both */ else if(minunstab==UNSTABLE) { p_stabp->bmax=p_stabp->bmin; p_stabp->bmin=p_stabp->bmin-p_stabp->increment; cut.b=p_stabp->bmin; sim.b=p_stabp->bmin; if (cut.bbminabs) { cut.b=p_stabp->bminabs; sim.b=p_stabp->bminabs; } init_param(&sim, &tool, &work, &cut,p_stabp); minunstab=p_schk(p_ypfun,p_rkcon,sim,p_hist,p_spres,&work, NO_PRINT,p_stabp); } /* If neither limit gave a good result change both until it works */ else if (minunstab==STAB_ERROR) { p_stabp->bmax=p_stabp->bmax+p_stabp->increment; cut.b=p_stabp->bmax; sim.b=p_stabp->bmax; if (cut.bbmaxabs) { cut.b=p_stabp->bmaxabs; sim.b=p_stabp->bmaxabs; } init_param(&sim, &tool, &work, &cut,p_stabp); maxunstab=p_schk(p_ypfun,p_rkcon,sim,p_hist,p_spres,&work, NO_PRINT,p_stabp); p_stabp->bmin=p_stabp->bmin-p_stabp->increment; cut.b=p_stabp->bmin; sim.b=p_stabp->bmin; if (cut.bbminabs) {

127 cut.b=p_stabp->bminabs; sim.b=p_stabp->bminabs; } init_param(&sim, &tool, &work, &cut,p_stabp); minunstab=p_schk(p_ypfun,p_rkcon,sim,p_hist,p_spres,&work, NO_PRINT,p_stabp); } } while ((p_stabp->bmax-p_stabp-> ((p_stabp->bmax-p_stabp->bmin)>p_stabp->tolerance) bmin)>p_stabp->tolerance) { bmean=(p_stabp->bmax+p_stabp->bmin)/2.0; cut.b=bmean; sim.b=bmean; printf("%f\n",bmean); init_param(&sim, &tool, &work, &cut,p_stabp); meanunstab=p_schk(p_ypfun,p_rkcon,sim,p_hist,p_spres,&work, NO_PRINT,p_stabp); if (meanunstab==0) return STAB_ERROR; else if (meanunstab*minunstab==-1) { p_stabp->bmax=bmean; maxunstab=meanunstab; } else { p_stabp->bmin=bmean; minunstab=meanunstab; } } return ((p_stabp->bmax+p_stabp->bmin ((p_stabp->bmax+p_stabp->bmin)/2.0); )/2.0); } int stab_check1(void (*p_ypfun)(),RK45CONST *p_rkcon,SIM_PARAM sim, HISTORY *p_hist,STEPDATA *p_spres,WORKPIECE_INFO *p_work, int printout,STAB_CHART_PARAM *p_stabp) /* This function checks the stability of a system by examining the growth or decay of the motions of the system. This function is used for for constant speed systems. */ { HISTORY *p_hfinal; int nstart,nfinal,changes,resul nstart,nfinal,changes,result; t; double xposdiff,xposdiff_last,period xposdiff,xposdiff_last,period; ; char *stepdatafile="stepdata.tx *stepdatafile="stepdata.txt"; t"; nstart=0; changes=0;

/* Number of consecutive same sign changes in amplitude */

/* Initialize the system by simulating INIT_ROTATIONS turns of the workpiece */ sim.tfinal=p_stabp->initrot*p_work->rotperiod; nstart=0; init_step(p_spres); result=rk45ad(p_ypfun,p_rkcon,&sim, p_hist,nstart, &p_hfinal,&nfinal, result=rk45ad(p_ypfun,p_rkcon,&sim, p_spres,NO_PRINT,NO_WRITEFILE,NO_APPEND,stepdatafile); /* Exits if cut is interrupted */ if (sim.int_cut!=0) return UNSTABLE;

128 xposdiff_last=0.0; /* Loops until rms vibration amplitude has either increased or decreased MAX_CHANGES revolutions of the workpiece */ while ((abs(changes)m ((abs(changes)maxchanges)&&(sim.totstepsrotperiod; /* Setup up the simulation to start where the last bit left off */ p_hist=p_hfinal; nstart=nfinal; /* Simulates the system for one revolution */ result=rk45ad(p_ypfun,p_rkcon,&sim,p_hist,nstart,&p_hfinal,&nfinal, p_spres,NO_PRINT,NO_WRITEFILE,APPEND,stepdatafile); /* Returns UNSTABLE if the cut was interrupted */ if (sim.int_cut!=0) return UNSTABLE; xposdiff=sim.xmax-sim.xmin; /* Determines some sort of mean value of the amplitude of the motion */ if (printout) printf("maxxpos= %e\n", sim.xmax); /* Compares the amplitude of the motion to that seen on the previous workpiece rotation */ if (xposdiff==0.0) return STABLE; if (xposdiff>xposdiff_last) { if (changes>=0.0) changes++; else changes=1; } else { if (changes<=0.0) changes--; else changes=-1; } xposdiff_last=xposdiff; } if (changes==p_stabp->maxchange (changes==p_stabp->maxchanges) s) return UNSTABLE; /* Returns one for unstable system */ else if (changes==-p_stabp->maxchang (changes==-p_stabp->maxchanges) es) return STABLE; /* Returns -1 for stable system */ else { printf("An error has occurred in the stability checking routine\n");

129 printf("Please contact your local Crapsoft (tm) representative \n"); printf("for details on how we would be glad to solve this problem\n"); printf("for a nominal fee.\n"); return STAB_ERROR; /* Returns 0 for a function error */ } } int stab_check2(void (*p_ypfun)(),RK45CONST *p_rkcon,SIM_PARAM sim, HISTORY *p_hist,STEPDATA *p_spres,WORKPIECE_INFO *p_work, int printout,STAB_CHART_PARAM *p_stabp) /* This function checks the stability of a system by examining the growth or decay of the motions of the system. This function is useful for systems systems in which the vibrations fall below a given threshold, but do not drop to zero */ { HISTORY *p_hfinal; int nstart,nfinal,changes,resul nstart,nfinal,changes,result,growing,stop; t,growing,stop; double xposdiff,xposdiff_last; char *stepdatafile="stepdata.tx *stepdatafile="stepdata.txt"; t"; changes=0; /* Number of consecutive same sign changes in amplitude */ stop=0; /* Initialize the system by simulating INIT_ROTATIONS turns of the workpiece */ sim.tfinal=p_stabp->initrot*p_work->rotperiod; nstart=0; init_step(p_spres); result=rk45ad(p_ypfun,p_rkcon,&sim, p_hist,nstart, &p_hfinal,&nfinal, result=rk45ad(p_ypfun,p_rkcon,&sim, p_spres,NO_PRINT,NO_WRITEFILE,NO_APPEND,stepdatafile); /* Exits if cut is interrupted */ if (sim.int_cut!=0) return UNSTABLE; xposdiff_last=0.0; growing=0; /* Loops until vibration amplitude has remained in the allowable band MAX_CHANGES revolutions of the workpiece */ while (!(stop)&&(sim.totstepsrotperiod; /* Simulates the system for one revolution */ result=rk45ad(p_ypfun,p_rkcon,&sim,p_hist,nstart,&p_hfinal,&nfinal, p_spres,NO_PRINT,NO_WRITEFILE,APPEND,stepdatafile); xposdiff=sim.xmax-sim.xmin;

130 /* Returns UNSTABLE if the cut was interrupted */ if (sim.int_cut!=0) return UNSTABLE; /* Determines some sort of mean value of the amplitude of the motion */ if (printout) printf("xposdiff= %e\n", xposdiff); /* Compares the amplitude of the motion to that seen on the previous workpiece rotation */ if (xposdiff==0.0) return STABLE; if (xposdiff>p_stabp->maxxmot (xposdiff>p_stabp->maxxmotion) ion) { if (changes>=0) changes++; else changes=1; if (changes==p_stabp->maxchanges (changes==p_stabp->maxchanges) ) stop=1; } else { if (changes<=0) changes--; else changes=-1; if (xposdiff>xposdiff_last) growing=growing+1; else growing=0; if (((-changes==p_stabp->maxchanges)&&(growingmaxc (((-changes==p_stabp->maxchanges)&&(growingmaxchanges))|| hanges))|| (growing>p_stabp->maxchanges)) { stop=1; if (growing>p_stabp->maxchange (growing>p_stabp->maxchanges) s) changes=p_stabp->maxchanges; } } xposdiff_last=xposdiff; } if (changes==p_stabp->maxchange (changes==p_stabp->maxchanges) s) return UNSTABLE; /* Returns one for unstable system */ else if (changes==-p_stabp->maxchang (changes==-p_stabp->maxchanges) es) return STABLE; /* Returns -1 for stable system */ else { printf("An error has occurred in the stability checking routine\n"); printf("Please contact your local Crapsoft (tm) representative \n"); printf("for details on how we would be glad to solve this problem\n"); printf("for a nominal fee.\n"); return STAB_ERROR; /* Returns 0 for a function error */ } } int stab_check3(void (*p_ypfun)(),RK45CONST *p_rkcon,SIM_PARAM sim, HISTORY *p_hist,STEPDATA *p_spres,WORKPIECE_INFO *p_work, int printout,STAB_CHART_PARAM *p_stabp) /* This function checks the stability of a system by examining the growth

131 or decay of the motions of the system. This function is useful for systems with variable spindle speed and constant feed rate */ { HISTORY *p_hfinal; int nstart,nfinal,result,i; double period; char *stepdatafile="stepdata.tx *stepdatafile="stepdata.txt"; t"; /* Sets up the simulation period in the event of variable speed machining */ period=1.0/sim.rpmfreq; sim.tfinal=period; nstart=0; init_step(p_spres); result=rk45ad(p_ypfun,p_rkcon,&sim, p_hist,nstart, &p_hfinal,&nfinal, result=rk45ad(p_ypfun,p_rkcon,&sim, p_spres,NO_PRINT,NO_WRITEFILE,NO_APPEND,stepdatafile); /* Exits if cut is interrupted */ if (sim.int_cut!=0) return UNSTABLE; /* Checks to see if amplitude stays within specified limits for MAX_CHANGES revolutions of the workpiece */ for (i=1;imaxchanges; (i=1;imaxchanges;i++) i++) { /* Initialize Stuff */ sim.xmax=-1000.0; sim.xmin=1000.0; sim.xvmax=-1000.0; sim.xvmin=1000.0; sim.xamax=-1000.0; sim.xamin=1000.0; /* Setup up the simulation to start where the last bit left off */ p_hist=p_hfinal; nstart=nfinal; /* Sets the final time so that the workpiece will spin only one time */ sim.tfinal=sim.tfinal+period; /* Simulates the system for one revolution */ result=rk45ad(p_ypfun,p_rkcon,&sim,p_hist,nstart,&p_hfinal,&nfinal, p_spres,NO_PRINT,NO_WRITEFILE,APPEND,stepdatafile); /* Returns UNSTABLE if the cut was interrupted */ if (sim.int_cut!=0) return UNSTABLE; /* Compares the amplitude of the motion to that seen on the previous workpiece rotation */ if (sim.xmax==sim.xmin) return STABLE; if ((sim.xmax>sim.xmaxdisp)|| ((sim.xmax>sim.xmaxdisp)||(sim.xmin
132 printf("for a nominal fee.\n"); return STAB_ERROR; /* Returns 0 for a function error */ } } return STABLE; } int stab_check4(void (*p_ypfun)(),RK45CONST *p_rkcon,SIM_PARAM sim, HISTORY *p_hist,STEPDATA *p_spres,WORKPIECE_INFO *p_work, int printout,STAB_CHART_PARAM *p_stabp) /* This function checks the stability of a system by examining the growth or decay of the motions of the system. This function is useful for systems with variable spindle speed and constant feed rate */ { HISTORY *p_hfinal; int nstart,nfinal,result,i; double period,xaccdiff,xaccdifflast; char *stepdatafile="stepdata.tx *stepdatafile="stepdata.txt"; t"; /* Sets up the simulation period in the event of variable speed machining */ period=1.0/sim.rpmfreq; sim.tfinal=period; nstart=0; init_step(p_spres); result=rk45ad(p_ypfun,p_rkcon,&sim, p_hist,nstart, &p_hfinal,&nfinal, result=rk45ad(p_ypfun,p_rkcon,&sim, p_spres,NO_PRINT,NO_WRITEFILE,NO_APPEND,stepdatafile); xaccdifflast=sim.xamax-sim.xamin; /* Exits if cut is interrupted */ if (sim.int_cut!=0) return UNSTABLE; /* Checks to see if amplitude decays for MAX_CHANGES revolutions of the workpiece */ for (i=1;imaxchanges; (i=1;imaxchanges;i++) i++) { /* Initialize Stuff */ sim.xmax=-1000.0; sim.xmin=1000.0; sim.xvmax=-1000.0; sim.xvmin=1000.0; sim.xamax=-1000.0; sim.xamin=1000.0; /* Setup up the simulation to start where the last bit left off */ p_hist=p_hfinal; nstart=nfinal; /* Sets the final time so that the workpiece will spin only one time */ sim.tfinal=sim.tfinal+period; /* Simulates the system for one revolution */ result=rk45ad(p_ypfun,p_rkcon,&sim,p_hist,nstart,&p_hfinal,&nfinal, p_spres,NO_PRINT,NO_WRITEFILE,APPEND,stepdatafile); /* Returns UNSTABLE if the cut was interrupted */ if (sim.int_cut!=0) return UNSTABLE;

133 xaccdiff=sim.xamax-sim.xamin; /* Compares the amplitude of the motion to that seen on the previous workpiece rotation */ if (sim.xmax==sim.xmin) return STABLE; if (xaccdiff>xaccdifflast) { return UNSTABLE; } if (sim.totsteps==sim.maxiter (sim.totsteps==sim.maxiter) ) { printf("An error has occurred in the stability checking routine\n"); printf("Please contact your local Crapsoft (tm) representative \n"); printf("for details on how we would be glad to solve this problem\n"); printf("for a nominal fee.\n"); return STAB_ERROR; /* Returns 0 for a function error */ } } return STABLE; }

A.13. envelope.c /* envelope.c Thomas O'Reilly October 1, 1997 Fuinctions to record maximum values */ #include #include #include #include #include #include #include

"orthsim.h" "cutlist.h" "stabchk.h" "integ.h"

/* Function Prototypes */ int envelope(void (*p_ypfun)(),RK45CONST *p_rkcon,SIM_PARAM sim, HISTORY *p_hist,STEPDATA *p_spres, WORKPIECE_INFO *p_work, FILE *fp); void envsim(void (*p_ypfun)(),SIM_PARAM sim,STEPDATA *p_spres, HISTORY *p_hist,TOOL_INFO tool,WORKPIECE_INFO work, CUT_PARAM cut); /* Real stuff starts here */ void envsim(void (*p_ypfun)(),SIM_PARAM sim,STEPDATA *p_spres, HISTORY *p_hist,TOOL_INFO tool,WORKPIECE_INFO work, CUT_PARAM cut) { /* Variable type declaration */ int nfinal, result; HISTORY p_hfinal; FILE *fp; RK45CONST *p_rkcon; /* Structure of rk45 constants */

134 /* Allocates space for rk45 constants */ p_rkcon=(RK45CONST *) malloc(sizeof(RK45CONST)); /* Initializes rk45 constants */ rk45init(p_rkcon); fp=fopen("envelope.txt","w"); result=envelope(*p_ypfun,p_rkcon,sim,p_hist,p_spres,&work,fp); fclose(fp); if (result==-1) { printf("The simulation printf("but didn't run result=sim.maxiter; } else if (result>0) printf("The cut became else printf("The simulation

reached the maximum number of time steps\n"); to the desired number of rotations.\n");

interrupted during rotation %d\n",result); was successful.\n");

return; } int envelope(void (*p_ypfun)(),RK45CONST *p_rkcon,SIM_PARAM sim, HISTORY *p_hist,STEPDATA *p_spres,WORKPIECE_INFO *p_work, FILE *fp) /* This function checks the stability of a system by examining the growth or decay of the motions of the system */ { HISTORY *p_hfinal; int nstart,nfinal,result,i,rota nstart,nfinal,result,i,rotations; tions; char *stepdatafile="envelope.tx *stepdatafile="envelope.txt"; t"; rotations=floor(sim.tfinal/p_work->rotperiod); nstart=0; sim.tfinal=0.0; init_step(p_spres); fprintf(fp,"0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n"); /* Loops until rms vibration amplitude has either increased or decreased MAX_CHANGES revolutions of the workpiece */ for (i=0;i<=rotations;i++) { sim.xmax=-1000.0; sim.xmin=1000.0; sim.xvmax=-1000.0; sim.xvmin=1000.0; sim.xamax=-1000.0; sim.xamin=1000.0; sim.ymax=-1000.0; sim.ymin=1000.0; sim.yvmax=-1000.0; sim.yvmin=1000.0; sim.yamax=-1000.0; sim.yamin=1000.0; /* Sets the final time so that the workpiece will spin only one time */ sim.tfinal=sim.tfinal+p_work->rotperiod;

135 /* Simulates the system for one revolution */ result=rk45ad(p_ypfun,p_rkcon,&sim,p_hist,nstart,&p_hfinal,&nfinal, p_spres,NO_PRINT,NO_WRITEFILE,APPEND,"stepdata.txt"); /* Setup up the simulation to start where the last bit left off */ nstart=nfinal; p_hist=p_hfinal; fprintf(fp, "%e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e\n", sim.txmax,sim.xmax,sim.txmin,sim.xmin, sim.txvmax,sim.xvmax,sim.txvmin,sim.xvmin, sim.txamax,sim.xamax,sim.txamin,sim.xamin, sim.tymax,sim.ymax,sim.tymin,sim.ymin, sim.tyvmax,sim.yvmax,sim.tyvmin,sim.yvmin, sim.tyamax,sim.yamax,sim.tyamin,sim.yamin); printf("Totstep= %d \n", sim.totsteps); /* Compares the amplitude of the motion to that seen on the previous workpiece rotation */

// //

if (!result) { return -1; } if (sim.int_cut!=0) return i;

} return 0; }

136

Appendix B

SIMULATION CONFIGURATION

B.1. Configuration File mode=1 /* Program mode kx=5500000.0 /* Tool stiffness in x direction, N/m*/ ky=3300000.0 /* Tool stiffness in y direction, N/m*/ mx=0.48 /* Effective mass in x direction, kg*/ my=0.7 /* Effective mass in y direction, kg*/ zetax=0.009 /* Damping ratio for x direction, percent*/ zetay=0.0013 /* Damping ratio for y direction, percent*/ cutstiff=1e9 /* Workpiece cutting stiffness N/m^2 */ d0=0.050 /* Workpiece Diameter m */ N=600.0 /* Rotational Speed, RPM */ DOC=0.00001 /* Depth of cut, meters per revolution, m */ b=0.0005 /* Width of cut, m */ db_dt=0.0 /* Rate of change of width of cut , m/s*/ theta=20.0 /* Angle of cutting force , degrees */ epsmin=0.00000005 /* lower bound on error before increasing epsmax=0.0000001 /* Maximum error per step */ maxiter=2000000 /* Maximum number of time steps */ fileskip=5 /* Print time step data every fileskip steps */ tfinal=6.0 /* End time, s */ hinc=1.25 /* Ratio by which step size is increased if hdec=0.75 /* Ratio by which step size is decreased if hmaxrat=0.5 /* Ratio of maximum step size to initial hminrat=2.0 /* Ratio of minimum step size to initial bmax=0.01 /* Maximum starting width of cut, m bmaxabs=0.01 /* Maximum final width of cut, m bmin=0.0 /* Minimum starting width of cut, m bminabs=0.0 /* Minimum final width of cut, m slots=32.0 /* increment to shift upper and lower limit by tolerance=0.0000025 /* Tolerance on width of cut, m Nmax=2000.0 /* Maximum spindle speed, RPM Nmin=250.0 /* Minimum spindle speed, RPM Nsteps=100 /* Number of increments in spindle speed initrot=10 /* initial rotations for stability checking maxchanges=5 /* Number of changes required for stability test maxxmotion=1e-6 /* Maximum amplitude of x direction motion, m rpmfreq=0.0 /* Frequency of spindle speed changes, Hz rpmamp=0.0 /* Amplitude of spindle speed changes, % of mean speed constfeed=1 /* Whether the feed rate is constant or varies

137 B.2. Configuration Procedure

The file in the previous section is a sample of the configuration file used to set system and simulation parameters. This is a plain text file file named simulate.cfg that should be present in the directory from which the program is to be run. run. At the start of   program execution this file is read and the values it contains are used in place of the default values set in the executable. If the file is not found, the program will write write one using those default values. This file can then be edited as needed. Each line in the file defines a different parameter. Each line begins with a variable name that defines the variable to be set. This is followed by and equals sign sign and then the value. Each value is followed by at least one space, and then a comment describing the variable. There must be no spaces prior to the value, and there must be a space between the value and any comments. If a variable name is repeated the program will use the second value in the file.

B.3. Configuration Variables

Although most of the values in the file are adequately explained by their comment the following require more detailed explainations. mode: mode: Sets simulation modes: Mode 1: Time domain simulation. Mode 2: Constant speed stability chart. Mode 4: Variable speed stability chart. Mode 6: Time domain simulation producing an envelope file. hinc/hdec: hinc/hdec: If the error estimate estimate on a time step exceeds the allowable error that time step is repeated with a smaller step size. The new step size is hdec times the the old value. If the the error estimate is smaller than the allowable step size the step size for the next step will be hinc time the previous value. hmaxrat/hminrat : The simulation guesses an initial time step that is 1/10 the period of  the highest frequency mode. The maximum step size is then the initial value multiplied multiplied  by hmaxrat , the minimum is the initial divided by hminrat .  slots:  slots: When producing a stability chart the difference between bmax and bmin is divided  by slots  by slots.. If at some point during the process the stability limit limit at a given speed is not  bracketed by the upper and lower limit in use at the time both limits are raised or lowered  by this quotient until one is bracketed.

138 rpmamp: rpmamp: For variable speed simulations simulations the mean spindle speed, N  speed, N , is multiplied by this quantity to determine the amplitude of spindle speed oscillation. constfeed : If this equals 1, the feed feed rate will be constant, if this equals zero, the feed rate will vary with spindle speed.  Nsteps:  Nsteps: When producing a stability chart, the number of spindle speeds between Nmin between Nmin and Nmax and Nmax to find stability at.

139

Appendix C

SIMULATION OUTPUT FILE FORMATS

C.1. stepdata.txt

This file is produced or time domain simulations using mode 1. It is a plain text file with 10 columns. Column 1: Time, s Column 2: x 2: x-direction -direction displacement, m Column 3: x 3: x-direction -direction velocity, m/s Column 4: y 4: y-direction -direction displacement, m Column 5: y 5: y-direction -direction velocity, m/s 2 Column 6: x 6: x-direction -direction acceleration, m/s 2 Column 7: y 7: y-direction -direction acceleration, m/s Column 8: Chip thickness, m Column 9: Distance traveled by the cutting tool along the workpiece surface, m Column 10: Distance traveled by the cutting tool into the workpiece surface, m Column 9 and 10 are used primarily to debug the program.

C.2. stabch.txt

This file contains the results results of any stability chart generation (mode (mode 2 or 4.) It is a  plain text file with 2 columns. Column 1: Spindle speed, RPM Column 2: Stability limit, m

140 C.3. envelope.txt

This file contains the results results of time domain simulation simulation using mode 6. It is a plain text file with 24 columns. These columns contain the the upper and lower chatter envelopes for position, velocity and acceleration in the x the x and y and y directions and the time at which they occurred. Column 1: Time vector for column 2 Column 2: Upper envelope, x envelope, x-direction -direction displacement, m Column 3: Time vector for column 4 Column 4: Lower envelope, x envelope, x-direction -direction displacement, m Column 5: Time vector for column 6 Column 6: Upper envelope, x envelope, x-direction -direction velocity, m/s Column 7: Time vector for column 8 Column 8: Lower envelope, x envelope, x-direction -direction velocity, m/s Column 9: Time vector for column 10 2 Column 10: Upper envelope, x envelope, x-direction -direction acceleration, m/s Column 11: Time vector for column 12 2 Column 12: Lower envelope, x envelope, x-direction -direction acceleration, m/s Column 13: Time vector for column 14 Column 14: Upper envelope, y envelope, y-direction -direction displacement, m Column 15: Time vector for column 16 Column 16: Lower envelope, y envelope, y-direction -direction displacement, m Column 17: Time vector for column 18 Column 18: Upper envelope, y envelope, y-direction -direction velocity, m/s Column 19: Time vector for column 20 Column 20: Lower envelope, y envelope, y-direction -direction velocity, m/s Column 21: Time vector for column 22 2 Column 22: Upper envelope, y envelope, y-direction -direction acceleration, m/s Column 23: Time vector for column 24 2 Column 24: Lower envelope, y envelope, y-direction -direction acceleration, m/s

141

REFERENCES

1.

de Canniere, J. H. van Brussel and J. van Bogaert. Contribution to the Mathematical Analysis of Variable Speed Machining. Applied Mathematical Modeling, Vol. 5, No. 3, June 1981, pp158-164.

2.

Elbestawi, M.A., F. Ismail, R. Du and B.C. Ullagaddi. Modeling Machining Dynamics Including Damping in the Tool-Workpiece Too l-Workpiece Interface, Journal of  Engineering for Industry, vol. 116, Nov 1994, pp. 435-439.

3.

Hanna, N. H. and S. A. Tobias. A theory of Nonlinear Regenerative Chatter. Journal of Engineering for Industry, ASME Transactions, February 1974, pp 247255.

4.

Hsu, C. S. Application of the Tau-Decomposition Method to Dynamical Systems Subjected to Retarded Follower Forces. Journal of Applied Mechanics, June 1970, pp 259-266.

5.

Inamura, T. and T. Sata. Stability Analysis of Cutting under Varying Spindle Speed. CIRP Annals, vol. 23, 1974, pp. 119-120.

6.

Kato, S, E., E. Marui, M. Hashimoto and T. Yamada. Characteristics of Chatter  Vibration of Lathe Tools. Memoirs of the Faculty of Engineering, Nagoya University, v. 38, n. 2, 1986, pp. 208-215.

7.

Kato, S, E., E. Marui, M. Hashimoto and T. Yamada. Characteristics of  Dynamics Cutting Force and Mechanism of Chatter Vibration. Memoirs of the Faculty of Engineering, Nagoya University, v. 38, n. 2, 1986, pp. 216-222.

8.

Kato, S, E., E. Marui, M. Hashimoto and T. Yamada. Analytical Considerations of Chatter Vibration of Lathe Tools. Memoirs of the Faculty of Engineering,  Nagoya University, v. 38, n. 2, 1986, pp. 223-232.

9.

Lee, M. S., and C. S. Hsu. On the τ-Decomposition Method of Stability Analysis for Retarded Dynamical Systems, SIAM Journal on Control, Vo l. 7, 1969, pp 242-259.

10.

Ota, Hiroshi, Eiji Kondo and Toshikatsu Yamada. Regenerative Re generative Chatter  Vibrations of Turning Workpieces. Workpieces. JSME International Journal, Ser. III, III, Vol. 32,  No. 1, 1989, pp 142-149.

142 11.

Rhorer, Richard L. and Don A. Lucca. Applocations of Cutting Mechanics Principles to Ultraprecision Machining and Grinding, Vol. 1: Session Notes. Tutorial notes, 1997 APSE Annual Meeting, Norfolk, VA, October 1997.

12.

Sexton, J. S. and B. J. Stone. Investigation of the Transient Effects During Variable Speed Cutting. Cutting. Journal of Mechanical Engineering Science, Vol. 22,  No. 3, June 1980, pp 107-118.

13.

Sexton, J. S. and B. J. Stone. Stability of Machining with Continuously Varying Spindle Speed. Annals of CIRP, Vol. 28, No. 1, 1978, pp 321-326.

14.

Sexton, J. S. R. D. Milne and B. J. Stone. Stability Analysis of Single Point Machining with Varying Spindle Speed. Applied Mathematical Modeling, Vol. 1, No. 6, Sept. 1977, pp310-318.

15.

Smith, S. T. and D. G. Chetwynd. Foundations of Ultraprecision Mechanism Design. Gordon and Breach Science Publishers, 1992.

16.

Stepan, G. Retarded dynamical systems: stability and characteristic functions, Longman Scientific & Technical , Essex, UK, 1989.

17.

Stone, B. J. Effect of Time Varying Parameters on the Build-Up of Chatter in Turning. Annals of CIRP, Vol. 34, No. 1, 1985, pp 371-374.

18.

Takemura, T. and T. Kitamura, T. Hoshi. Active Suppression of Chatter by Programmed Variation of Spindle Speed. CIRP Annals, v 23, 1974, pp. 121-122.

19.

Tlusty, J. Machine Dynamics (1985). In King, Robert I. (ed), Handbook of High Speed Machining Technology. (pp. (pp. 48-143). New York, Chapman and Hall.

20.

Tsao, Tsu-Chin, Mark W. McCarthy and Shiv G. Kapoor. A New Approach to Stability Analysis of Variable Speed Machining Systems. International Journal of  Machine Tools and Manufacturing. vol. 33, no. 6, pp. 791-808, 1993.

21.

Yantek, David. Chatter and Process Damping in Diamond Turning. M. S. Thesis, The Pennsylvania State University, December 1996.