Modeling and Control of Cutter Vibration in the Presence of Random Distribution of Microhardness of Workpiece and Nonlinear Cutting Forces in Lathe Process

This work presents a comprehensive approach to the control of tool's position, in the presence of machine tool structure vibration, nonlinear cutting force, and random tool vibration due to random distribution of microhardness of workpiece material. The controller is combination of Proportional and linear quadratic gaussian- (P-LQG-) type constructed from an augmented model of both tool-actuator dynamics and a nonlinear dynamic model relating tool displacement to cutting forces. The latter model is obtained using black-box system identification of experimental orthogonal cutting data in which tool displacement is the input and cutting force is the output. The controller is evaluated and its performance is demonstrated.


Introduction
Performance of machining processes is assessed by dimensional and geometrical accuracy as well as the surface texture of the part. Factors such as cutting force, workpiece vibration, machine-tool vibration, process instability, tool wear, and thermal deformation deteriorate this performance. Recent interest in high-speed machining and high manufacturing efficiency requires faster, higher bandwidth actuation and controllers that are more robust. Besides determining part functional behavior, the surface texture also plays a key role in the area of manufacturing process control [1]. The basis for using surface texture as means of process control is derived from the fact that a slight change in the manufacturing process manifests itself as a corresponding change in the resulting surface geometry. The roughness, waviness, lay, and flaws constitute the texture of the work piece surface [2]. Surface texture of the work piece is highly dependent on tool position during cutting which in turn is affected by cutting dynamics. Therefore, extensive research work has been done on modeling and analysis of tool workpiece interaction and the dynamics of the cutting process. Altintas et al. [3] have presented a cutting force model as a function of regenerative chip thickness, cutting speed, and the velocity and acceleration terms of the vibrations. In their work, amplitude and frequency of inner and outer vibration waves are generated on the chip by an instrumented fast tool servo-powered by a piezo actuator. They established a model in which the coefficients of dynamic cutting forces are identified and used in analyzing the effect of cutting speed, tool wear, vibration frequency, and wavelength on the chatter stability of a turning process. Ikua et al. [4] have studied the cutting force generated in ball end milling and presented a model for the radial and tangential forces as a function of depth of cut and other factors related to the geometry of the tool and workpieces.
Researchers have tried to establish quantitative relations between the surface finish and the cutting parameters of the depth of cut, feed, and cutting speed using methods such as wavelength decomposition surface roughness, wavelet analysis, and tool vibration [4][5][6][7][8][9]. The cutterworkpiece interaction forces are assumed to cause relative displacement between the cutter and the work piece, which influences the surface generation mechanism. Analytical, experimental and mechanistic methods were used to predict the interaction forces. The analytical methods focused on establishing a relationship between the cutting force and the instantaneous uncut chip cross-section [10,11] and the nonlinear mechanisms [12,13]. The analytical models were not capable of predicting the dynamic forces accurately due to the secondary nonlinear effects that stems from the tool/workpiece interaction. The experimental methods include static, dynamic, and time-series methods. The static methods are assumed based on linear assumptions and linear cutting conditions [14]. The dynamic methods replaced the static methods, and the cutting process is assumed to be a combination of two independent actions, namely, the wave cutting action and the wave removing action [15]. The time-series method was formulated to identify the dynamic cutting force coefficients as well as the transfer functions of the three-dimensional dynamic cutting process directly from operating data [16]. The mechanistic modeling methods view the machining process as a combination of the chip load/cutting force relationship, cutting-tool geometry, machining conditions, and tool/work-piece displacement due to cutting forces. The workpiece dynamic behavior is constructed using either the distributed-parameter [17] or the lumped-parameter approaches [18]. Few researchers investigated the dynamic effect of the tool holder. Shawky and Elbestawi [19] concluded that the effect of tangential force is insignificant, and no model of the tool holder is needed in that direction. Their experimental work showed that modeling the tool by a second order system in each axis provides significant approximation of the dynamic behavior of a single-point tool during cutting [19].
Control of machining process includes three levels of control that one might encounter in a controller for machining process; namely, they are the servo control, the process control [20], and the supervisory control [21]. In the servo control process, the motion of the cutting tool is taken relative to that of the workpiece. The process control level is used to control the cutting forces and tool wear to maintain high production rates and good part quality. The highest level of control is the supervisory control and it directly measures product-related variables such as part dimensions and surface finish. Different approaches are used to correct for machining errors by means of dynamic error compensation. More recently, Tian et al. [22] have presented a methodology for modeling and control of a high-precision flexure-based mechanism for ultraprecision turning operation based on the position control of an auxiliary precision mechanism utilized on the turret of the conventional lathe to implement nanometer level in feed. El-Sinawi [23] has presented an optimal control method for controlling the tool position in orthogonal cutting in both feed are radial directions. Moradi et al. [24] have presented a robust control method for orthogonal turning process in the process was modeled as a single degree of freedom model that includes the effect of tool flank wear with a control input of the system being force applied to the tool provided by a piezoactuator. Huang et al. [25] have presented a tool wear detection based on cutting force monitoring.
In this work, a control scheme for the purpose of improving surface texture of turned surfaces through control of tool position in the presence of nonlinear cutting force is developed.
Improvement of surface quality is achieved via active positioning of the tool or cutter through the attenuation of cutting forces effect on the tool in both radial and tangential directions. The process is assumed to be stochastic due to both process and measurement noise. Experimental force-displacement data is used to construct a nonlinear dynamic model of dynamic cutting forces. Various nonlinear models are constructed using system identification techniques including ARX, nonlinear ARX, and Hammerstein-Wiener techniques [26]. The best nonlinear model obtained that closely fits experimental data is then linearized and constructed in state-space form to later utilize in the linear quadratic Gaussian controller. Simulation results will be used to verify the effectiveness of the proposed modeling and control approach and to enhance current understanding of the mechanisms responsible for generating both stochastic and deterministic components of the surface texture. Simulation work of the proposed modeling and control techniques will be based on actual parameters for actuators and the cutting process data, which is readily available in the literature [27].

Nonhomogeneous Distribution of Microhardness.
A statistical approach was proposed by Zhang and Kapoor [28] to formulate the nonhomogeneous distribution of microhardness which is based on random sampling processes in mixing of chemicals presented by Tucker [29]. A representative part is chosen from the work-piece material. The microstructural constituents as well as their microhardness values could be experimentally identified, say H i for i = 1, 2, . . . , n, where n represents the number of measurement locations on the chosen part. Based on these data, the three statistical parameters of μ a , variance σ 2 a , and the correlation coefficient function ρ(r), are calculated using the standard statistical formulas. Note that the parameter r of the correlation coefficient function represents the distance in space between the ith and i th measured locations. These three parameters provide quantitative information on the size, the shape, and the segregation of microstructures and depict the microhardness distribution present in the material being machined. They also form a basis for quantifying the effect of microconstituent scale on the bulk properties of hardness of the material being machined. Figure 1, the geometric shape of each sample is a parallelepiped being πD/n s long and a cross-sectional area equal to a product of depth of cut and feed. Due to the nonhomogeneity, each sample holds its own distinct ratio of microstructures. Therefore, μ s is different from sample to sample. The sample hardness means μ s j for the jth sample can be obtained as

Sampling Process and Parameter Estimation. As shown in
where N s is the number of locations within a sample (a subset of measurement locations, n) and j = 1, 2, . . . , n s . The mean hardness value of all sample (μ s1 , μ s2 , . . . , μ sns ) would be equal if the microstructures in the material were uniformly distributed (in a microscale). However, this is never the case in practice [28]. As a result, the sample hardness mean, μ s , differs from sample to sample. In fact, this variation reveals, from a microbase analysis, the main source of random excitation phenomenon observed during machining.
When the cutting tool is cutting a sample (a portion of the circumference of the work-piece being cut) with a small value of mean hardness, the magnitude of the cutting force drops. On the other hand, the cutting force jumps up when the cutting tool is meeting a sample with a large value of mean hardness. Because the variation in the mean hardness value from sample to sample is of random nature, the dynamic variation of the cutting force is also random in nature. The larger the mean hardness value variation, the stronger the random excitation system present during machining. Since the sample hardness mean μ s , the key parameter of the random excitation system, is a random variable, it is best dealt with statistically. It can be shown that the mean of the sample mean distribution is equal to the mean of its population distribution that is, μ s = μ a. . However, estimating the variance of the sample hardness mean distribution or the sample mean variance, using the standard statistical formulas, is not practical, because such an evaluation would require the subgrouping of the identified locations of microstructures into individual samples. Therefore, this procedure to evaluate σ 2 s would be extremely tedious if applied in practice. Furthermore, the main drawback of this procedure is that quantitative relationship between the sample mean variance and the geometric sample shape; namely, the three cutting parameters cannot be easily established. For this reason, the sample mean variance theory [29] is employed to establish this quantitative relationship. The sample mean variance theory stemmed from the need to characterize the composition uniformity of a chemical mixing process, a similar process of taking random samples from the population distribution and evaluating the sample mean values afterwards. The theory offers a mathematical approach to quantitatively evaluate the sample mean variance based on the knowledge of the mean, the variance, and the correlation coefficient function of the population distribution, from which the random samples are being taken.
To estimate the variance of the sample hardness mean distribution σ 2 s , we start by assuming the work-piece material to be a two-component mixture consisting of components A and B. This indicates that the cutting edge (tool) experience different microstructural portions during cutting. The analysis may be generalized into multiple-component mixtures by letting one component be A and all others (including free space) to be B. in this case the analysis must be performed (n − 1) times for a mixture of n components. If samples with identical geometric shape (each has a volume V ) are taken from the mixture, the average concentration of component A in an identical sample is given by while the average of μ s will depend on the average of a; that is, For convenience, we define the quantity c as so that the deviation of μ s from the average for any particular sample is The square deviation of μ s is now the double integral is taken over all pairs of incremental volumes within the sample. If the square deviation of μ s is averaged over all samples in the mixture, one gets the sample variance σ 2 Let M be the number of samples in the mixture. If we take y i to be the spatial location of a reference point in the ith sample, and if x is referred to a local coordinate system with its origin at y i and oriented with respect to the sample, then (6) and (7) can be combined to give The order of the integration and summation can be reversed, bringing the summation inside the integral. Then, inspection of (5) shows that where σ 2 a is the variance in concentration among all samples in the mixture, ρ(r) is called the correlation coefficient function, and r is the magnitude of the distance between the points x and x for Equation (8) can now be expressed as the average of the correlation coefficient function over a single representative sample The outer integral in (11) specifies an incremental volume dv, for which the inner integral is evaluated. This inner integral can be expressed in spherical coordinates, with the limits on r, θ, and φ chosen to limit the integral to points within the sample volume as shown in Figure 1. Note that r is simply the distance between the volume increment dv and the increment dv specified by the outer integral. Let x be the location of this later incremental volume. To avoid the difficulty of specifying limits on θ and φ in (12), we define a function W * (x, r) as follows. Consider a spherical shell of radius r and thickness dr, centered at x. W * (x, r) is the fraction of the volume of that shell lying inside the sample as shown in Figure 3. Since the integrand in (12), we define a function of r only, one can easily evaluate the integrals with respect to θ and φ when the integrand is multiplied by W * . That is, The infinite upper limit is used on r, because W * will be zero when r is so great that the entire shell lies outside the sample.
Substituting back into (11) and rearranging gives The terms within the square brackets in (14) are related to the sample geometry only. Define this part as the sample shape function W(x, r), because it can be integrated with respect to the sample volume separately.
by substituting (15) into (14), the formula to evaluate the sample mean variance is given by Substituting for V from Figure 1, (16) becomes Note that (17) consists of two parts. One part is the integral, and the other is the constant term (4n s σ 2 a /D f d), where n s is the number of samples, D is the workpiece diameter, f is the feed rate, and d is the depth of cut. For a given work-piece material under different cutting parameter settings, the constant term decreases as any of the terms in the denominator increases. The integral term always increases as any of the three parameters f, d, and/or D increases, because the geometric shape factor W(r) always increases (sometimes very slightly), and at the same time, the correlation coefficient function ρ(r) remains unchanged as long as the material of the work-piece remains unchanged. Therefore, the magnitude of the sample mean variance σ 2 s is a compromise between these two terms. In fact, (17) serves as a mathematical model to predict values of the sample mean variance under different cutting parameter settings.
It is well known in statistics that the sampling distributions of sample means are normal even though the samples may be drawn from a nonnormal population such as uniform distributions. The mathematical justification of this relationship is the central limit theorem. Based on this theorem, the distribution of the sample hardness mean, μ s , is normal. The mean of the sample hardness means is expressed as where μ a is the mean of the population distribution of hardness and the variance of the sample hardness means is given by (17).
Since the mean level of the dynamic cutting force is directly related to the mean value μ s , this indicates that the random excitation system (dynamic cutting force) can be mathematically describes by a normal distribution as shown in Figure 2.
The variance of the normal distribution σ 2 s represents the variation level o the dynamic cutting force about its mean level. The larger the variance σ 2 s , the more significant the μ s variation as well as the variation level of the dynamic cutting force. As the variance σ 2 s decreases, the variation level of the dynamic cutting force decreases accordingly. When the variance σ 2 s approaches zero, indicating that all sample mean hardness values are equal to each other, the random excitations caused by the cutting force diminishes.
Finally, the following very important conclusions can be drawn for the preceding sections.
(1) The random excitation system can be mathematically described by a normal distribution. (2) The two parameters μ s and σ 2 s of the normal distribution fully describe the characteristics of the random excitation system. (3) The variation level of the dynamic cutting force represented by σ 2 s , is directly affected by the variation of the cutting parameters; see (17).
These conclusions are vital to the design of the proposed Kalman estimator, which will be the presented later on.

Nonlinear Force Model
Force-displacement data obtained from experimental cutting of a 6061 Aluminum workpiece using a carbide tool with 0 • rake angle and 7 • clearance angle, depth of cut of 0.4 mm, feed of 0.050 mm/rev, and spindle speed: 2200 rev/min, which are shown in Figures 5 and 6, have been used to construct a force displacement model using black-box system identification.
Four dynamic models with tool displacement as input and cutting force as output are constructed are indentified and compared as shown in Figure 3. In reference to Figure 3, the first model identified as LinMod1 is a linear ARX model of the form where f c , x c are the cutting force and tool displacement in the corresponding direction, respectively.
are delayed input and output variables called regressors. Linear ARX model predicts the output f c as a weighted sum of its regressors. τ is the number of past output terms, while μ is the number of past input terms used to predict the current output f c , and e(t) is a white noise sequence. The second and third models identified as LinMod3 and LinMod4 are constructed as follows. LinMod3 is constructed using predictive error minimization or PEM with various constraints on the order or number of states of the state-space model. This method yields a discrete time-domain state-space dynamic model of the form (20) in which x are the states, y is output, A, B, C, K, D represent dynamics, input output, disturbance and feed through matrices, respectively [30]. Narx1 model is constructed using Nonlinear ARX model which is an extension of ARX model given in (1) except that the input-output mapping is nonlinear and of the form; where f n is the output and F is a nonlinear mapping function, such as wavelet network, tree partition, and so forth, [31]. Nhw1 model shown in Figure 3 is constructed using Hammerstein-Wiener model; see [31] for further details. It is clear from Figure 3 that the Narx1 has the best fit to experimental data with 99.9% fit.

Linearization of the Nonlinear ARX Model.
The nonlinear model that best fits input-output data is clearly the nonlinear ARX model mapped with a single-layer sigmoid function. The nonlinear mapping of input-output data done with where v is a vector of the regressors. L(v − r) + d is the output of the linear block, d is a scalar offset, g(Q(v−r)) is the output of the nonlinear function block, and Q is a projection matrix and r is the mean of regressors (v). Function g(v) expressed as g(v) = n k=1 α k φ(β k (v − ξ k )) with φ(s) = (e s + 1) −1 is the nonlinearity estimator of the sigmoid function type. β k is a Linearized state-space model of the radial force iṡ X r (t) = A r X r (t) + B r x cr , whereẊ r (t) are the states of the radial force model, A r , B r , C r , and D r are dynamic, input, output, and feed through matrices, respectively. f x , x cr are radial cutting force and radial tool displacement, respectively. The same procedure is carried out for the tangential cutting force model which yieldsẊ Figure 4 shows the difference between tool post (passive) and proposed active tool fixture that will allow for implementation of the proposed controller. Figure 4(a) shows the conventional tool post rigidly attached to the machine tool structure, while Figure 4(b) shows the active tool holder platform equipped with two actuators placed between the tool and the machine tool structure. This will allow the control force to provide necessary manipulation of the tool to maintain a constant depth of cut. During cutting, the tool is perturbed form nominal depth of cut by two inputs; namely, the dynamic cutting  forces in both radial and tangential directions. w(t) represents nonhomogenous micro-hardness distribution model presented in Section 2, producing the stochastic component of the cutting force.

Modeling of the Machining Process.
State-space model of the tool-actuator assembly in the radial direction can be represented aṡ x r = A cr x r + B cr F ar , Y cr = C cr x r + D cr F ar . (25) such that A cr = 0 1 −kac/m −bac/m , B cr = 0 1/m , C cr = [1 0], D cr = 0, where x r is vector of the tool states, namely, displacement and velocity. m is the mass of the tool and its mount, k ac , b ac are the elastic stiffness and damping coefficient of the actuator, respectively. k ac , b ac can be obtained by an FR test of the tool-actuator assembly. Notice that the same procedure can be carried out for the tangential direction, since both actuators are orthogonal. Refer to [23,27] for more information on the model development.

Control Strategy and the Controller Design
The control strategy will be centered on reducing the amplitude of the tool' dynamic displacement to zero for the purpose of maintaining a constant depth of cut and subsequently, a smooth surface texture of the workpiece. This requires minimization of the error between the desired and the actual position of the tool. The desired or nominal tool position is the one yielding a constant depth of cut. However, the existence of the dynamic cutting forces perturbs the tool form its nominal position, and thus, varies the depth of cut causing subsequent deterioration in surface texture. Attenuating the effect of the cutting force on the tool's position, the error in the tool position will be minimized. This can be achieved by reacting on the tool with an equal but opposite forces through proper actuation. This task is difficult due to (a) the existence of process and measurement noise and (b) the nonlinearity of the dynamic cutting forces. To overcome these difficulties, the controller must be able to estimate the force needed to minimize the error in the tool position from noisy process using measurement data contaminated with noise. In addition to that, the controller has to track and manipulate the position of the tool effectively in the presence of nonlinear dynamic cutting forces. Added to all of the above, the controller has to maintain high stability and performance under various disturbance characteristics. To do so, the controller has to be constructed based on a model that takes into account the dynamics of the tool-actuator system as well as the cutting force dynamics. Therefore, The plant used for this purpose is a combination (i.e., augmentation) of the dynamic cutting forces and toolactuator dynamics in both radial and tangential directions. The first step in designing the proposed active controller is the design of the LQG estimator and regulator gains L e and K R , respectively. The analysis presented here is in state-space such that the augmented system depicted by (23) and (25) is; where the matrices A a , B a , C a , and D a are, respectively, the dynamic, input, output, and direct transmission matrices of the augmented system. The four matrices in the foregoing are used to design the LQG estimator and regulator gain matrices L e and K R such that This augmentation is necessary to incorporate the cutting force dynamics in the process of determining the optimal estimator and regulator gains. Equation (27) shows that L e and K R are both partitioned in two parts each, where [] cr corrects x cr (i.e., tool-actuator model states) and [] r corrects X r (i.e., cutting force states). The LQG (virtual model of the system) should be subjected to all inputs that the actual plant is subjected to, including the control force. Figure 7 shows the control scheme implementation in the radial direction. The proportional controller is used to adjust for any variation in the cutting force due to variation of any or all of the cutting parameters.

Numerical Example
Force-displacement results obtained from experimental cutting of a 6061 Aluminum workpiece using a carbide tool with 0 • rake angle and 7 • clearance angle, depth of cut of 0.4 mm, feed of 0.050 mm/rev, and spindle speed: 2200 rev/min, shown in Figures 5 and 6, are used to construct a forcedisplacement model using black-box system identification. Linearized State-space model obtained by nonlinear ARX modeling presented here in canonical form iṡ x cr , ETREMA actuator with peak-to-peak excursion of 50 × 10 −6 m is utilized as the active manipulator of the tool. The actuator's elastic stiffness and damping used in this study were taken as follows [27], m = 0.53 kg for the moving toolactuator assembly. k ac = 14.6 × 10 6 N/m, and b ac = 10 kg/s. Figure 8 shows the base force transmitted from the machinetool structure to the cutter through the actuator. This force is due to the vibration of the machine tool at low resonance frequencies as indicated in [27]. Figure 9 shows the vibration of the cutter due to random distribution of microhardness according to the model developed in Section 2 of this paper and represented by w(t) in Figure 7. The measurement noise v(t) is assumed to be random with zero mean and variation of ±10% of measurement values.

Results and Analysis
State-space model of the augmented system given in (26) is used to determine the LQG gains. LQG gains are highly dependent on weight matrices of the performance index. See [32] for further details on the selection of performance index weights.
The controller is implemented in both radial and tangential directions as outlined in Figure 7. However, results here are only for the radial direction (depth of cut direction), since the implementation in the feed direction is only a duplicate of the radial. Figure 10 shows the tool's displacement in the radial direction with and without control, where it is clear from the latter that the controller has succeeded in reducing the dynamic tool-displacement by approximately 40%-50%. The controller maintained excellent stability in the presence of process and measurement noise as shown by the bounded response of the cutter after control implementation. Control force in the radial direction is shown in Figure 11, which is well within the actuator's force capabilities. Cutting force as modeled by the nonlinear system identification shows that the cutting force has a significant magnitude especially when the random force due to nonhomogenous distribution of micro-hardness is taken into consideration. Figure 12 shows the cutting force to be reduced significantly when the tool displacement (i.e., vibration) is reduced as shown in Figure 10. The proportional gain used in simulation is P = 1 × 10 6 and can be tuned for better performance of the control scheme with upper limit determined based on the saturation force of the actuator.

Conclusions
In this research work, active control of an orthogonal cutting process is presented. Nonlinear dynamic cutting force model is generated from actual machining data obtained in both radial and tangential directions. An LQG-based controller is developed based on an augmented model of both toolactuator dynamics and a linearized model of the dynamic cutting force. The control objective is to eliminate the tool's dynamic displacement and maintain a constant depth of cut.
Simulation results have shown that the proposed control strategy has managed to significantly reduce the dynamic displacement of the tool with minimal force and calculations efforts even in the presence of significant randomness in the process and measurement as well as the nonlinearity of the cutting force. Experimental study is undergoing to verify the integrity and performance of the controller in real-time applications.
are then subjected to
(A. 4) In view of the assumption of small disturbances, the cross terms with disturbance. The terms w(t) and v(t) in (A.3) above are white disturbances with the following covariance properties: , u * (t)), R 12 (t) = r(x * (t), u * (t))Ew(t)v T (t)m T (x * (t), u * (t)). (A.5) The model is now linear in the vicinity of nominal trajectory.