Computer numerical control vertical machining centre feed drive modelling using the transmission line technique

Abstract

This study presents a novel application of the Transmission Line Matrix Method (TLM) for the modelling of the dynamic behaviour of non-linear hybrid systems for computer numerical control (CNC) machine tool drives. The application of the TLM technique implies the dividing of the ball-screw shaft into a number of identical elements in order to achieve the synchronisation of events in the simulation, and to provide an acceptable resolution according to the maximum frequency of interest. This entails the use of a high performance computing system with due consideration to the small time steps being applied in the simulation. Generally, the analysis of torsion and axial dynamic effects on a shaft implies the development of independent simulated models. This study presents a new procedure for the modelling of a ball-screw shaft by the synchronisation of the axial and torsion dynamics into the same model. The model parameters were obtained with equipments such as laser interferometer, ball bar, electronic levels, signal acquisition systems, etc. The MTLM models for single and two-axis configurations have been simulated and matches well with the measured responses of machines. The new modelling approach designated the Modified Transmission Line Method (MTLM) extends the TLM approach retaining all its inherent qualities but gives improved convergence and processing speeds. Further work since, not the subject of this paper, have identified its potential for real-time application.

Background

Ford ¹ proposed that the design process of any high precision computer numerical control (CNC) machine tool must embrace error avoidance, error measurement, and error compensation techniques in order to achieve an ever increasing demand for greater accuracy and therefore a more stringent performance specification. The measurement of the repeatable time and spatial errors can allow compensation methods to be applied to the machine for correction of those errors provided sufficient resolution has been allowed for in the design, with the absolute limit on accuracy for a particular machine being defined by its measured repeatability figures and the responsiveness of its compensating axes. The main areas of concern affecting component accuracy are environmental effects, user effects and the machine tool static and dynamic accuracy. The sources of error confined to a machine tool are geometrical inaccuracies, thermally induced errors and load errors.

Machine tool feed drive modelling and system identification

Pislaru et al. ² stated that the modelling used for designing high performance CNC machine tool feed drives should accurately represent the system consisting of three cascaded closed loops for acceleration, velocity and position control. They utilised a modular approach and a hybrid modelling technique which included for distributed load, explicit damping coefficients, backlash and friction. The model was a combination of distributed and lumped elements described by partial differential equations and ordinary differential equations as proposed by Bartlett and Whalley. ³

Holroyd et al. ⁴ modelled the dynamic behaviour of a ball-screw with moving nut using a finite element approach. The ball-screw was divided into a large number of elements and contact/boundary conditions for each element and adjacent ones were studied. An important conclusion was that the natural frequencies of the ball-screw system vary in time due to two causes: (a) the lateral restraint produced by the nut when the screw was transversally vibrated and (b) the relation between screw torsion, axial motion and worktable/saddle tilt.

Transmission line matrix method

The transmission line matrix (TLM) method belongs to the general class of differential time domain numerical modelling methods. It is used to solve time-dependent or transient problems, thus involving ordinary and partial differential equations (PDE).

The application of TLM to the modelling of fluid systems by Boucher and Kitsios ⁵ and mechanical systems by Partridge et al. ⁶ implies the same representation used in the hybrid approach by Pislaru, ² i.e. lumped parameter modules, components with localised dynamic effects and distributed parameter modules. The difference resides on the fact that transmission line equations introduce natural time delays between components due to the wave propagation speed in a transmission line. Therefore, distributed components must be divided into a number of identical elements in order to achieve the synchronisation of events in the simulation, and to produce acceptable resolution according to the maximum frequency of the simulation. This characteristic gives the possibility to include the effect of the nut movement suggested in Holroyd’s et al. ⁴ approach.

Johns and Beurle ⁷ described the two equivalent circuits used in the TLM technique: the stub and the link circuits. TLM links are two-port one-dimensional building blocks that can be used for one-, two- or three-dimensional modelling. TLM stubs are one-port units which can be used for solving circuits and equations, and are used in multi-dimensional modelling to complement TLM links. Links are used for modelling distributed parameter elements and stubs for lumped parameter elements.

Christopoulis ⁸ stated that any electrical circuit could be represented as a network of transmission line circuits by simply replacing the reactive elements with corresponding stubs. Variables such as voltage and current are regarded as discrete pulses bouncing at a velocity to and from nodes of these stubs at each time step. The voltage and current in each element are determined from the incident and reflected pulses in each stub. He demonstrated the transmission line and its Thevenin equivalent.

The MTLM stub, taking into account changes of sample time was applied to the modelling of a rotor shaft arrangement by Partridge et al., ⁶ where Castaneda ⁹ demonstrated improved robustness and accuracy of the simulation. A reduction of 40% in the number of mathematical operations and almost 35% reduction in the mean square modelling error figures were achieved. This method of approach was utilised for this application.

Single axis CNC test rig

The single-axis test rig (Figure 1a) represents the y-axis of a vertical machining centre (VMC). The single axis is fitted with a Heidenhain TNC-426PB motion controller, a SIMODRIVE-611 Siemens inverter, and a ball-screw arrangement directly coupled to a permanent magnet synchronous motor.

Figure 1.

(a) Test rig and (b) test rig block model.

The Heidenhain ¹⁰ TNC-426PB motion controller offers digital control of functions such as interpolation, position control, speed control, current control and provides PWM signals to the inverter unit. The control of each axis is implemented with an algorithm following the principle of cascade control in conjunction with velocity and acceleration feed forward.

The SIMODRIVE-611 consists of a common feed module that provides the DC voltage link from the power supply mains and a set of drive modules that activate each motor. In the case of the test rig, the drive module consists of a power module (inverter) and an interface card that communicates with the TNC-426PB motion controller.

The nut of the ball screw system is preloaded and the screw shaft is mounted on preloaded bearings in a fixed–fixed configuration.

The model of the axis feed drive is built by interconnecting the TLM models of the rig elements. The dynamics acting on the axis feed drive are defined by the interrelation of three blocks: The motion controller, the inverter and motor electrical equations, and motor mechanical equations and mechanical transmission elements as illustrated in Figure 1(b). Table 1 contains the test rig specification.

Table 1.

Test rig specification.

Motion controller	HEIDENHAIN TNC 426PB (280 476-24)
Electrical drive	SIEMENS Simodrive 611 (6SN1123-1AA00-0CA1)
Motor	SIEMENS 1FT6082-1AF7-1AG1
Guide ways	THK SNS45-LC
Front bearing	Double Row Angular Contact Ball Bearing ZKLF 3080.2RS PE
Rear bearing	Double Row Angular Contact Ball Bearing ZKLF 3080.2RS PE
Ballscrew	THK BIF (BNFN) 4016-5 RRG0S-1346(1065)L-C3-E
Coupling	ROBA – ES 28 940.000
Rotary encoder	HEIDENHAIN ERN 1387
Linear encoder	HEIDENHAIN AE LS486C ML620
Transverse	500 mm
Rapid Traverses	30,000 mm/min
Feed Rates	100 – 12000 mm/min
Table Size	115 mm × 58 mm
Table mass	312 kg
Saddle mass	524 kg
Load mass	854 kg

TLM method applied to CNC machine tool feed drive

Motion controller

Heidenhain ¹⁰ operates with position and velocity control in cascade with feed forward capability as shown in Figures 2 and 3. Seborg et al. ¹¹ concluded that feed forward control improves trajectory tracking compensating for the effect of disturbances before they affect the response of the controlled system. The velocity feed forward path connects the reference position to the velocity loop through the gain, k_ff (Figure 3). When the position profile changes, the velocity feed forward transfers immediately that change to the velocity command. This speeds up the system response, thus not relying solely on the position loop. The primary shortcoming of velocity feed forward is that it induces overshoot when changes on the velocity direction occur. Acceleration feed forward eliminates the overshoot caused by velocity feed forward without reducing loop gains. The basic principle of feed forward control is to inject the position and velocity set points and their correspondent first derivatives at appropriate points in the control loops. Heinemann and Papiernik ¹² demonstrated that the following errors are minimised, thus achieving a significant increase in the dynamic response to position and velocity set point changes.

Figure 2.

Motion controller block model.

Figure 3.

Motion controller flow model.

Altinas ¹³ demonstrated that the interpolator (Figure 2) calculates a velocity profile v (t), according to defined motion commands: feed rate (f_r), maximum acceleration (a_max), maximum jerk (j_max), and desired displacement (d). The velocity profile is then used to generate the set of reference positions in order to achieve a smooth velocity transition.

Koren ¹⁴ and Bullock ¹⁵ demonstrated that linear and circular interpolation, respectively, are performed to keep the tool path velocity constant at the given feed rate along a straight line in a plane of motion.

If the system has multiple axes, the motions of individual axes are co-ordinated. When the individual axes are not co-ordinated, each of the axes will start moving at the same time, but finish at separate times producing slew motion as defined by Hugh. ¹⁶

The interpolator overcomes this effect using an interpolation technique. Position values from the interpolator can be filtered in order to obtain a smoother motion profile. The filtered values constitute the reference position signals for the position controllers.

The Heidenhain TNC426PB ¹⁷ motion controller (Figure 3) was implemented in software featuring the algorithms for position, velocity and current control at different sampling rates. ¹⁰ That is:

The interpolator generates a reference position value d_e every 3 ms.

The position controller generates a reference velocity value v_ref every 3 ms.

The velocity controller generates a reference current value i_qref every 0.6 ms.

The current controller gives a reference voltage value e_dqref to the PWM generator at a rate of 0.2 ms.

The dynamic response-matching filter (first order delay filter) was used to delay the position profile signal according to the transient response during acceleration and deceleration (the equivalent position time constant of the closed position control loop). Delay values can be set in the interval of 1–255 ms.

The jerk limitation filter was used to adapt the position profile to the machine dynamics in order to attain high machining velocity. The coefficients of the filter were calculated according to the minimum order of the filter and the tolerance for contour transitions.

The PT₂, second order element was used to include a delay in the reference current (i_ref) to damp the frequency interference oscillations (circa: 0.3–2 ms).

Sliding friction was compensated within the range of the velocity controller by compensating the sliding friction at low velocity and at the rated velocity of the motor. The compensation at low velocity was achieved by feeding forward the reference current value (measured at approximately 10 rpm) at every change in direction. The compensation at the rated velocity was achieved by feeding forward the current i_{ff_rv} according to the value of the reference velocity (equation (2)). A delay filter was included to prevent overcompensation when the traverse direction was reversed at high feed rates.

$i_{\mathrm{ff}}(k)=i_{\mathrm{ff}\_l}+i_{\mathrm{ff}\_\mathrm{rv}}(k)$

(1)

$i_{\mathrm{ff}\_\mathrm{rv}}(k)=\{i_{\mathrm{ffm}}{v}_{\mathrm{ref}}(k)\ge i_{\mathrm{ffm}}{v}_{\mathrm{ref}}(k)k_{\mathrm{iff}}\mathrm{for}\hspace{1em}|v_{\mathrm{ref}}(k)|<i_{\mathrm{ffm}}-i_{\mathrm{ffm}}\hspace{1em}{v}_{\mathrm{ref}}(k)\le -i_{\mathrm{ffm}}$

(2) where

• i_ffm is the reference current measured at the rated velocity of the motor; i_{ff_l} is the reference current at 10 rpm and k_iff is the scaling factor.

Position controller with velocity feed forward

The position controller evaluates the difference between the reference and actual position value d_act (from the rotary or linear encoder) to calculate the position error d_e. A reference velocity value (v_ref) is then generated every t_s seconds according to the following

$v_{\mathrm{ref}}=d_e{k}_v+v_{\mathrm{ff}}$

(3)

$d_e=d_{\mathrm{ref}}-d_{\mathrm{act}}$

(4)

$v_{\mathrm{ff}}=k_{\mathrm{rff}}\frac{d}{\mathrm{dt}}{d}_{\mathrm{ref}}$

(5) The signal d_act is calculated by equation (6) if the linear encoder is used as a position feedback system. Equation (7) gives the value for d_act when the rotary encoder is used instead

$d_{\mathrm{act}}(k)=d_l(k)$

(6)

$d_{\mathrm{act}}(k)=k_b{\theta }_m(k)$

(7) where

• k_b is the force to torque constant of the ball screw.

The TLM model for this module is represented by equations (3) to (7) for a sample time t_s = 3 ms.

By a combination of Christopoulis ⁸ and Castaneda ⁹ theory

$v_{\mathrm{ff}}(k)=d_{\mathrm{ref}}(k)Z_{\mathrm{vff}}+E_{\mathrm{vff}}^i(k)$

(8)

$E_{\mathrm{vff}}^i(k+1)=-d_{\mathrm{ref}}(k)Z_{\mathrm{vff}}$

(9)

$v_{\mathrm{ref}}(k)=k_v{d}_e(k)+v_{\mathrm{ff}}(k)$

(10)

$d_e(k)=d_{\mathrm{ref}}(k)-d_{\mathrm{act}}(k)$

(11) where $Z_{\mathrm{vff}}=k_{\mathrm{vff}}/t_s$ and Z_vff is the characteristic impedance of the stub and $E_{\mathrm{vff}}^i$ (k) is the incident pulse.

Velocity controller with acceleration feed forward

The velocity controller evaluates the difference between the reference and actual velocity value (v_act) measured with a rotary encoder attached to the motor. A reference current value (i_ref) is generated every t_s s by a Proportional-Integral-Differential (PID) strategy.

Acceleration feed-forward is used in parallel with the velocity controller in order to minimise the spikes caused by changes in velocity direction. A signal i_ht (holding current) is added to the acceleration feed forward to counterbalance the axis load when the axis in the vertical position.

Equations (13) to (18) represent the TLM model for this module (t_s= 0.6 ms). The following relationship calculates the reference current signal

$i_{\mathrm{ref}}=i_p+i_{\text{int}}+i_{\mathrm{der}}+a_{\mathrm{ff}}+i_{\mathrm{ht}}$

(12) where

$i_p=k_p{v}_e$

(13)

$i_{\mathrm{int}}=k_i\int v_e\text{d}t$

(14)

$i_{\mathrm{der}}=k_d\frac{{\mathrm{dv}}_e}{\mathrm{dt}}$

(15)

$a_{\mathrm{ff}}=k_{\mathrm{aff}}\frac{{\mathrm{dv}}_{\mathrm{ff}}}{\mathrm{dt}}$

(16)

$v_e=v_{\mathrm{ref}}-v_{\mathrm{act}}$

(17) Here k_p, k_i and k_d are the proportional, integral and derivative gains of the velocity error v_e, respectively. k_aff is the acceleration feed forward gain. Applying the Castaneda ⁹ TLM transform (equations (19) to (26)) for the cycle time t_s gives

$i_{\mathrm{ref}}(k)=i_p(k)+i_{\text{int}}(k)+i_{\mathrm{der}}(k)+a_{\mathrm{ff}}(k)$

(18)

$i_p(k)=k_p{v}_e(k)$

(19)

$i_{\text{int}}(k)=Z_i{v}_e(k)+E_i^i(k)$

(20)

$E_i^i(k+1)=i_{\text{int}}(k)$

(21)

$i_{\mathrm{der}}(k)=Z_d{v}_e(k)+E_d^i(k)$

(22)

$E_d^i(k+1)=-Z_d{v}_e(k)$

(23)

$a_{\mathrm{fft}}(k)=Z_{\mathrm{aff}}{v}_{\mathrm{ff}}\hspace{1em}(k)+E_{\mathrm{aff}}^i(k)$

(24)

$E_{\mathrm{aff}}^i(k+1)=-Z_{\mathrm{aff}}{v}_{\mathrm{ff}}\hspace{1em}(k)$

(25)

$v_e(k)=v_{\mathrm{ref}}(k)-v_{\mathrm{act}}(k)$

(26) where

$Z_i=k_i{t}_s;\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}{Z}_d=k_d/t_s\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}\text{and}\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}{Z}_{\mathrm{aff}}=k_{\mathrm{aff}}/t_s$

Band-stop filter

Heidenhain ¹⁷ states that a filter is generally used in the velocity feedback to damp the fundamental frequency of the control system, when it is higher than 500 Hz. A first-order low-pass filter is used when the oscillation frequency is between 500 and 700 Hz. A second-order low-pass filter is used if the oscillation frequency is higher than 700 Hz.

When the controlled system is insufficiently damped, it will be impossible to attain a sufficiently short settling time without inducing oscillations in the step response of the velocity controller. The step response will oscillate even with a low proportional factor. A second-order lag element (PT₂) is used to include a delay in the reference current (i_ref) to damp the frequency interference oscillations.

A band-rejection filter is included in series with the PT₂ element to damp oscillations that cannot be compensated by the differential factor of the velocity controller, the PT2 element or the low-pass filter.

Mathworks ¹⁸ state these digital filters are generally modelled as implementations of the standard difference equation

$a(1)i_{\mathrm{qref}}(k)\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}=b(1)i_{\mathrm{ref}}(k)+b(2)i_{\mathrm{ref}}(k-1)+\dots \mathrm{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}}+b(2N+1)i_{\mathrm{ref}}(k-2N)-a(2)i_{\mathrm{qref}}(k-1)-\dots \mathrm{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}}-a(2N+1)i_{\mathrm{qref}}(k-2N)$

(27) where

• b and a represent the numerator and denominator coefficients of the filter and N is the order of the filter.

The band-stop filter was implemented as indicated by Mathworks ¹⁹ in the transposed direct-form II structure of equation (27). This is a canonical form that has the minimum number of delay elements. At sample k, the routine computes the difference equations as follows

$i_{\mathrm{qref}}(k)=\mathrm{num}(1)i_{\mathrm{ref}}(k)+{\mathrm{zz}}_1(k-1)$

(28)

${\mathrm{zz}}_1(k)=\mathrm{num}(2)i_{\mathrm{ref}}(k)+{\mathrm{zz}}_2(k-1)-\mathrm{den}(2)i_{\mathrm{qref}}\hspace{1em}(k)$

(29)

${\mathrm{zz}}_{\mathrm{len}-3}(k)=\mathrm{num}(\mathrm{len}-2)i_{\mathrm{ref}}(k)+{\mathrm{zz}}_{\mathrm{len}-2}(k-1)\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}-\mathrm{den}(\mathrm{len}-2)i_{\mathrm{qref}}(k)$

(30)

${\mathrm{zz}}_{\mathrm{len}-2}(k)=\mathrm{num}(\mathrm{len}0)i_{\mathrm{ref}}(k)+{\mathrm{zz}}_{\mathrm{len}0}(k-1)\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}-\mathrm{den}(\mathrm{len}0)i_{\mathrm{qref}}(k)$

(31)

${\mathrm{zz}}_{\mathrm{len}0}(k)=\mathrm{num}(\mathrm{len})i_{\mathrm{ref}}(k)-\mathrm{den}(\mathrm{len})i_{\mathrm{qref}}(k)$

(32)

$\mathrm{len}=\mathrm{len}0+1$

(33) where len0 is the filter order, and num and den represent the numerator and denominator filter coefficients. The delay outputs zz_i(1), $i=1,\dots ,\mathrm{len}0$ are initialised to 0. This is equivalent to assuming that both past inputs and outputs are zero.

Current controller

The current controller is an implementation of the model when the integral term k_ci is equal to zero.

Castaneda ⁹ developed the model shown and converted to the TLM transform of equations which are presented as

$i_{\mathrm{de}}(k)=i_{\mathrm{dref}}(k)-i_{\mathrm{dact}}(k)$

(34)

$e_{\mathrm{dref}}(k)=k_{\mathrm{cp}}{i}_{\mathrm{de}}(k)$

(35)

$i_{\mathrm{qe}}(k)=i_{\mathrm{qref}}(k)-i_{\mathrm{qact}}(k)$

(36)

$e_{\mathrm{qref}}(k)=k_{\mathrm{cp}}{i}_{\mathrm{qe}}(k)$

(37)

Inverter and motor (electrical)

Moreton ²⁰ describes the voltage relationships and provided the general equation for a permanent magnet synchronous motor (PMSM). Simon ²¹ explained the Clarke /Park transformation and their inverse relating the vector currents used for the generation of electromagnetic torque (Figure 4a).

Figure 4.

TLM model for the motor and coupling. (a) Motor electrical equations, (b) TLM model of the motor and coupling, (c) reduced TLM model and (d) Thevenin equivalent circuit.

The equations of performance for a PMSM in the d–q coordinate system representation are

$L_d\frac{d}{\mathrm{dt}}{i}_d=e_d-{\mathrm{Ri}}_d+L_q{\omega }_e{i}_q$

(38)

$L_q\frac{d}{\mathrm{dt}}{i}_q=e_q-{\mathrm{Ri}}_q-L_d{\omega }_e{i}_d-\lambda {\omega }_e$

(39) in which λ signifies the amplitude of the flux induced by the permanent magnets of the rotor in the stator phases; L_d and L_q are the d and q axis inductances. Subscripts d and q are used to identify the direct and quadrature currents/voltages, respectively. The generated electromagnetic torque T_e is defined as

$T_e=1.5p(\lambda i_q+(L_d-L_q)i_d{i}_q)$

(40) if the current $i_d\to 0$ , $T_e\to 1.5p\lambda i_q$ ; hence, the generated electromagnetic torque can be approximated using the torque constant of the motor, k_T

$T_e=k_T{i}_q$

(41) Castaneda ⁹ considered a motor with a round rotor and a simulation time step t_s, and derived the TLM transforms of equations (38), (39) and (41) giving

$Z_L{i}_d(k)+E_d^i(k)=e_d(k)-{\mathrm{Ri}}_d(k)+L{\omega }_e(k)i_q(k)$

(42)

$\text{where}:Z_L=L/t_s$

(43)

$E_d^i(k+1)=-i_d(k)Z_L$

(44)

$Z_L{i}_q(k)+E_q^i(k)=e_q(k)-{\mathrm{Ri}}_q(k)-{\omega }_e(k)({\mathrm{Li}}_d(k)+\lambda )$

(45)

$E_q^i(k+1)=-i_q(k)Z_L$

(46)

$T_e(k)=k_T{i}_q(k)$

(47) Rearranging equations (43) and (45) provides

$e_d(k)-E_d^i(k)-(Z_L+R)i_d(k)=-{\omega }_e(k){\mathrm{Li}}_q(k)$

(48)

$e_q(k)-E_q^i(k)-(Z_L+R)i_q(k)={\omega }_e(k)({\mathrm{Li}}_d(k)+\lambda )$

(49) For simulation purposes, it can be assumed that the terms on the left side of equations (48) and (49) are the resultants of the operations on the right side of the equations. Thus, these equations are rewritten as

$i_d(k)=M_{\mathrm{dq}}(e_d(k)-E_d^i(k)+{\omega }_e(k-1){\mathrm{Li}}_q(k-1))$

(50)

$i_q(k)=M_{\mathrm{dq}}(e_q(k)-E_q^i(k)-{\omega }_e(k-1)({\mathrm{Li}}_d(k-1)+\lambda ))$

(51) where

$M_{\mathrm{dq}}=1/(Z_L+R)$

(52)

Motor mechanical equations and coupling

The inertias of the rotary encoder and the coupling (hub 1) can be added to the inertia of the motor to simplify the calculations. Hence

$J_{\mathrm{mc}}=J_m+J_c+J_{\mathrm{re}}$

(53) where

• J_mc = motor plus coupling inertia; J_m = motor inertia and J_re = rotary encoder inertia.

This reduction of the model is possible because those inertias are modelled as lumped parameter elements. The resultant TLM model is illustrated in Figure 4(c) (reduced from Figure 4b).This electric circuit is solved finding the Thevenin equivalent ⁸ with respect to T_cs (Figure 4d), thus

$E_{\mathrm{eq}}(k)=Z_{\mathrm{Ecs}}{E}_{\mathrm{cs}}^i(k)+Z_{\mathrm{Ect}}{E}_{\mathrm{ct}}(k)$

(54)

$Z_{\mathrm{eq}}=\frac{Z_{\mathrm{cs}}{Z}_{\mathrm{ct}}}{Z_{\mathrm{cs}}+Z_{\mathrm{ct}}}$

(55) where

$Z_{\mathrm{Ecs}}=Z_{\mathrm{ct}}/(Z_{\mathrm{cs}}+Z_{\mathrm{ct}})$

(56) and

$Z_{\mathrm{Ect}}=Z_{\mathrm{cs}}/(Z_{\mathrm{cs}}+Z_{\mathrm{ct}})$

(57) and

$Z_{\mathrm{ct}}=Z_c+Z_t$

(58)

$E_{\mathrm{ct}}(k)=E_c^i(k)+2A_1^i(k)$

(59)

${\omega }_m(k)=M_{\omega m}(T_e(k)-E_{\mathrm{mc}}^i(k)-E_{\mathrm{eq}}(k))$

(60)

$T_{\mathrm{cs}}(k)={\omega }_m(k)Z_{\mathrm{eq}}+E_{\mathrm{eq}}(k)$

(61)

${\omega }_1(k)=M_{\omega 1}(T_{\mathrm{cs}}(k)-E_{\mathrm{ct}}(k))$

(62)

$E_{\mathrm{cs}}^i(k+1)=T_{\mathrm{cs}}(k)$

(63)

$E_c^i(k+1)=-Z_c{\omega }_1(k)$

(64)

$E_{\mathrm{mc}}^i(k+1)=-Z_{\mathrm{mc}}{\omega }_m(k)$

(65)

$B_1^i(k+1)={\omega }_1(k)Z_t+A_1^i(k)$

(66) where

$M_{\omega m}=1/(b_m+Z_{\mathrm{mc}}+Z_{\mathrm{eq}})$

(67) and

$M_{\omega 1}=1/Z_{\mathrm{ct}}$

(68) and

$Z_{\mathrm{mc}}=J_{\mathrm{mc}}/t_t$

(69) Equations (54) to (69) represent the TLM model of the motor mechanical equations and coupling.

Screw shaft torsional model

From the theory outlined by Partridge et al., ⁶ Holroyd et al., ⁴ and Castaneda, ⁹ the following model was derived.

The presence of the supporting bearings in the TLM model of the shaft generates the reflection of pulses arriving to the sections where they are placed. In that case, the propagation of pulses in the TLM model takes place on two specific zones (loops) as it is graphically represented in Figure 5(a). The front bearing is placed on section f_b, the nut is on section n_t and the rear bearing is on section h_t, where

$h_t=\mathrm{round}(l_r/l_{\mathrm{tor}})$

(70)

$l_{\mathrm{end}}=l_{\mathrm{ss}}-h_t{l}_{\mathrm{tor}}$

(71)

$J_{\mathrm{end}}=l_{\mathrm{end}}{I}_o$

(72)

$f_b=\mathrm{round}(l_f/l_{\mathrm{tor}})$

(73)

$n_t=\mathrm{ceil}(l_n/l_{\mathrm{tor}})$

(74) The inclusion of the nut in the model will cause the reflection of pulses arriving to section n_t, and therefore splitting the zone 2 in two loops as shown in Figure 5(b). The model is then reduced to the calculation of the angular velocity on section one, f_b, n_t and h_t; and the propagation of pulses on the other sections.

Figure 5.

Pulses propagation model for the screw shaft torsional model with moving nut. (a) Without nut and (b) Zone 2 including the nut.

The velocity of the front bearing (ω_fb+1) is calculated including the TLM model derived for the bearing's friction ²² (see Figure 6(a)), thus

$T=2(B_{\mathrm{fb}}^i(k)-A_{\mathrm{fb}+1}^i(k))$

(75)

${\omega }_{\mathrm{fb}+1}(k)\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}=\{0(T-\mathrm{sign}(T)T_{\mathrm{fb}1})/Z_{\mathrm{fbeq}}\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}\mathrm{for}\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}|T|<T_{\mathrm{fb}1}|T|\ge T_{\mathrm{fb}1}$

(76) where

$Z_{\mathrm{fbeq}}=2Z_t+b_{\mathrm{fb}}$

(77)

$A_{\mathrm{fb}}^i(k+1)=B_{\mathrm{fb}}^i(k)-{\omega }_{\mathrm{fb}+1}(k)Z_t$

(78) next pulses

$B_{\mathrm{fb}+1}^i(k+1)=A_{\mathrm{fb}+1}^i(k)+{\omega }_{\mathrm{fb}+1}(k)Z_t$

(79) The angular velocity ω_nt + 1 is calculated according to the procedure specified in the screw shaft, nut and table section.

Figure 6.

(a) Section f_b of the torsional model. (b) Section h_t of the torsional model.

The pulse propagation is defined by

$B_{\mathrm{nt}+1}^i(k+1)=A_{\mathrm{nt}+1}^i(k)+{\omega }_{\mathrm{nt}+1}(k)Z_t$

(80)

$A_{\mathrm{nt}}^i(k+1)=B_{\mathrm{nt}}^i(k)-{\omega }_{\mathrm{nt}+1}(k)Z_t$

(81) The velocity of the rear bearing (ω_ht + 1) is calculated using the procedure applied to the front bearing (see Figure 6b).

$T=2B_{\mathrm{ht}}^i(k)-E_{\mathrm{Jend}}^i(k)$

(82)

${\omega }_{\mathrm{ht}+1}(k)\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}=\{0(T-\mathrm{sign}(T)T_{\mathrm{rb}1})/Z_{\mathrm{rbeq}}\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}\mathrm{for}\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}|T|<T_{\mathrm{rb}1}|T|\ge T_{\mathrm{rb}1}$

(83)

$Z_{\mathrm{rbeq}}=Z_t+b_{\mathrm{rb}}+Z_{\mathrm{Jend}}$

(84) next pulses

$A_{\mathrm{ht}}^i(k+1)=B_{\mathrm{ht}}^i(k)-{\omega }_{\mathrm{ht}+1}(k)Z_t$

(85)

$E_{\mathrm{Jend}}^i(k+1)=-{\omega }_{\mathrm{ht}+1}(k)Z_{\mathrm{Jend}}$

(86) The propagation of $A^i$ and $B^i$ pulses on the other sections are given by

$B_j^i(k+1)=B_{j-1}^i(k)\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}\text{for}j=2,\dots ,h_t\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}j\ne f_b+1,n_t+1$

(87)

$A_j^i(k+1)=A_{j+1}^i(k)\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}\text{for}\hspace{1em}j=1,\dots ,h_t-1\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}j\ne f_b,n_t,h_t$

(88)

Screw shaft axial model

From the theory outlined by Partridge et al. ⁶ and Holroyd et al., ⁴ Castaneda ⁹ derived the following model.

The presence of the supporting bearings and the nut generates the reflection of pulses arriving to the sections where they are placed. This leads to the propagation of pulses on one zone (loops) in the axial model as it is graphically represented in Figure 7.

Figure 7.

Pulses propagation model for the screw shaft axial model with moving nut. (a) Without nut and (b) including the nut.

The inclusion of the nut in the model will cause the reflection of pulses arriving to section n_a, and therefore splits the model in two loops as shown in Figure 7(b). The front bearing is placed on the first section, the nut is on section n_a and the rear bearing is on section h_a, where

$h=\mathrm{round}(l_r/l_{\mathrm{axial}})$

(89)

$l_{\mathrm{end}}=l_{\mathrm{ss}}-l_{\mathrm{axial}}h$

(90)

$f_{\mathrm{ba}}=\mathrm{round}(l_f/l_{\mathrm{axial}})$

(91)

$l_{\mathrm{front}}=f_{\mathrm{ba}}{l}_{\mathrm{axial}}$

(92)

$h_a=(h-f_{\mathrm{ba}})$

(93)

$n_a=\mathrm{ceil}((l_n-l_{\mathrm{front}})/l_{\mathrm{tor}})$

(94) The model is reduced to the calculation of the velocities v_1a, v _{na  + 1} and v _{ha  + 1}; and the pulse propagation on the other sections as defined by the procedure presented by Castaneda ⁹ showing the derivation of the velocity v_1a; velocity v _{ha  + 1} and the pulse propagation.

Screw shaft, nut and table

Figure 8 illustrates the connection of the axial and torsional TLM screw shaft models with the nut and table models, thus

${\omega }_{\mathrm{nt}+1}(k)=\{0M_{\mathrm{wnt}1}(T-\mathrm{sign}(T)T_p)\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}\mathrm{for}\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}|T|\le T_p|T|>T_p$

(95)

$v_{\mathrm{nt}+1}(k)=k_b{\omega }_{\mathrm{nt}+1}(k)$

(96)

$T_a(k)=k_b{F}_a(k-1)$

(97) where

$T=2\mathrm{BA}-T_a(k)$

(98)

$\mathrm{BA}=B_{\mathrm{nt}}^i(k)-A_{\mathrm{nt}+1}^i(k)$

(99)

$M_{\mathrm{wnt}1}(k)=1/(2Z_t)$

(100) The components of the frictional force F_gw are calculated as follows

$F_{\mathrm{Ry}}=m_{\mathrm{ly}}+F_{\mathrm{cz}}$

(101)

$F_{\mathrm{Ly}}=|F_{\mathrm{cx}}|$

(102)

$F_E=X_y{F}_{\mathrm{Ry}}+Y_y{F}_{\mathrm{Ly}}$

(103)

$F_0=F_{\mathrm{gw}0}+b_{\mathrm{gw}}{F}_E$

(104)

$F_1=b_{\mathrm{gw}}{v}_l$

(105) The velocity of the table (v_l) is calculated finding the Thevenin equivalent ⁸ with respect to F_a (Figure 8c):

$E_{\mathrm{eq}}=v_{n+1}(k)Z_{\mathrm{eq}}+Z_{\mathrm{CD}}\mathrm{CD}+Z_{\mathrm{Ens}}{E}_{\mathrm{ns}}^i(k)$

(106)

$\mathrm{CD}=B_{\mathrm{na}}^i(k)-A_{\mathrm{na}+1}^i(k)$

(107) where

$Z_{\mathrm{eq}}=2Z_a{Z}_{\mathrm{ns}}/(2Z_a+Z_{\mathrm{ns}})$

(108)

$Z_{\mathrm{CD}}=2Z_{\mathrm{ns}}/(2Z_a+Z_{\mathrm{ns}})$

(109)

$Z_{\mathrm{Ens}}=2Z_a/(2Z_a+Z_{\mathrm{ns}})$

(110)

$v_l(k)=\{0M_{\mathrm{vl}}(F-\mathrm{sign}(F)F_0)\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}\mathrm{for}\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}|F|\le F_0|F|>F_0$

(111)

$F=E_{\mathrm{eq}}-E_l^i(k)-F_{\mathrm{cy}}(k)$

(112) where

$M_{\mathrm{vl}}=1/(Z_{\mathrm{eq}}+b_{\mathrm{gw}}+Z_l)$

(113)

$F_a(k)=E_{\mathrm{eq}}-v_l(k)Z_{\mathrm{eq}}$

(114)

$E_l^i(k+1)=-v_l(k)Z_l$

(115)

$E_{\mathrm{ns}}^i(k+1)=F_a(k)$

(116)

$v_{\mathrm{na}+1}(k)=M_{\mathrm{vna}1}(2\mathrm{CD}-F_a(k))$

(117)

$M_{\mathrm{vna}1}=1/(2Z_a)$

(118)

$d_l(k)=v_l(k)t_a+E_{\mathrm{dl}}^i(k)$

(119)

$E_{\mathrm{dl}}^i(k+1)=d_l(k)$

(120) It must be noted that a new ballscrew with preload is assumed to have no or minimal backlash. Preload is lost over time due to wear of balls and ball races and non-inclusion of backlash term would reduce the value of the model. The model is therefore made applicable to both cases:

• ball-screw with pretension in the nut (Figure 8) Backlash = 0 and

• ball-screw without pretension in the nut (Figure 9) T_p = 0 and Backlash ≠ 0.

The backlash model presented has been reduced to the following two possible states:

when the screw shaft is not in contact with the nut (Figure 9a) and

when the screw shaft is in contact with the nut (Figure 9b).

The state in which the axis will start at the beginning of the simulation depends on the following conditions:

• Non-contact if:

• –d_d = d_a and the direction of motion is negative (nut moving towards the motor).

• –d_d = d_a + Backlash and the direction of motion is positive.

Contact if:

–d_d = d_a and the direction of motion is positive.

–d_d = d_a + Backlash and the direction of motion is negative.

Variables T_a and F_a are zero if the screw shaft is not in contact with the nut. Thus the values for the velocities are given by the following equations

${\omega }_{\mathrm{nt}+1}(k)=\mathrm{BA}/Z_t$

(121)

$v_{\mathrm{nt}+1}(k)=k_b{\omega }_{\mathrm{nt}+1}(k)$

(122)

$v_{\mathrm{na}+1}(k)=\mathrm{CD}/Z_a$

(123)

$v_a(k)=v_{\mathrm{nt}+1}(k)+v_{\mathrm{na}+1}(k)$

(124)

$v_n(k)=-E_{\mathrm{ns}}^i(k)/Z_{\mathrm{ns}}$

(125)

$v_l(k)=\{0M_{\mathrm{vl}}(F-\mathrm{sign}(F)F_0)\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}\mathrm{for}\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}|F|\le F_0|F|>F_0$

(126)

$v_d(k)=v_l(k)+v_n(k)$

(127)

$F=-E_l^i(k)-F_{\mathrm{cy}}(k)$

(128) where

$M_{\mathrm{vl}}=1/(b_{\mathrm{gw}}+Z_l)$

(129)

$E_{\mathrm{ns}}^i(k+1)=0$

(130)

$F_a(k)=0$

(131)

$E_l^i(k+1)=-v_l(k)Z_l$

(132) The positions d_a and d_d are calculated by integrating the velocities for the sampling time t_a (equations (124) and (127)). d_l is calculated as in the case for preloaded nut. The model remains in this state if 0 < d_d(k) − d_a(k) < Backlash; otherwise the model changes to the contact state (Figure 9 with T_p = 0). The model will switch to the non-contact state when the sign of velocity v_nt + 1 changes.

Figure 9.

Backlash model. (a) Non-contact and (b) in contact.

Figure 8.

TLM model of the connection between nut and screw shaft. (a) Torsional model connextion, (b) axial load connextion and (c) Equivalent model for the connection between nut and screw shaft.

Implementation of two-axis TLM models to a CNC machine tool

Figure 10(c) shows the three-axis vertical machining centre where the unidirectional Z-axis maximum tool travel is 560 mm, and the bidirectional machining tables X-axis maximum travel and Y-axis maximum travel are both 510 mm. The Siemens 840D CNC controls all the machine tool functions. The spindle speed can be controlled from 60 RPM (revolutions per minute) to 8000 RPM and the maximum axis feed-rates are 12 m/min (fast traverse 30 m/min). The machine utilises Siemens 611 electrical feed drives with permanent magnet synchronous motors (Types 1FK6063 on X- and Y-axis and Type 1FT6064 on Z-axis). The motors for each axis are flexibly coupled to a ball-screw supported by bearings at both ends and utilise linear bearing guide-ways. The Siemens spindle drive unit provides four quadrant operations through a separately excited Siemens DC motor providing a constant torque range 60 RPM to 1500 RPM and above 1500 RPM constant power range (9 KW) to the maximum speed of 8000 RPM.

Figure 10.

Machine tool. (a) x–y axis feed drive, (b) block diagram and (c) machine tool and coordinate system.

The VMC-500 is fitted with a Siemens SINUMERIK 840D motion controller, which commands the CNC kernel functions for interpolation and position control. The motion controller is connected to the drives and I/O units via a PROFIBUS-DP interface as shown by Siemens. ²³ Each axis integrates a SIMODRIVE-611 Siemens inverter and a ball-screw arrangement directly coupled to a permanent magnet synchronous motor. The ball-screw systems incorporate a preloaded nut and the screw shaft mounted on a fixed-supported bearing configuration.

The SIMODRIVE-611 unit consists of a common feed module that provides the DC voltage link from the power supply mains and a set of drive modules that activate each motor. Every drive module consists of a power module (inverter) and a closed-loop plug-in unit. The closed-loop plug-in unit is dedicated to velocity control, current control and PWM generation functions. Table 2 contains the technical data for the VMC-500 Machine centre.

Table 2.

VMC-500 machine specifications.

Specification	X-axis	Y-axis	Z-axis
Motion Controller	SIEMENS SINUMERIK 840D (S7-315-2DP)
Electrical Drive	SIEMENS SIMODRIVE 611D (6SN1123-1AA00-0CA1)
Motor	SIEMENS 1FK6063-6AF71-1AG2		SIEMENS 1FT6064-6AH71-4AB1
Guideways	THK HSR35 A2 SS CO QZ + 1090 L H II
Front bearing	Double Row Angular Contact Ball Bearing NSK BSB 025062
Rear bearing	Single-row deep groove ball bearing NSK 6305 TB	Single-row radial ball bearing NSK 6005 – 2RS	Single-row deep groove ball bearing NSK 6305 TB
Ballscrew	THK
Coupling	ROTEX 24 GS spider 98 Shore A
Rotary encoder	HEIDENHAIN ERN 1387
Linear encoder	HEIDENHAIN AE LS186C ML620
Transverse	510 mm
Rapid Traverses	30000 mm/min
Feed Rates	300–12,000 mm/min
Table Size	700 mm × 520 mm
Table mass	327 Kg
Saddle mass		160 Kg
Carrier mass			212 Kg
Motor mass			152 Kg
Load mass			364 Kg

The main differences between the TNC 426PB and the SINUMERIK 840 D ²⁴ in terms of modelling its behaviour are as follows:

The SINUMERIK 840D includes a velocity response matching filter (first-order delay-filter) used to delay the velocity feed forward signal according to the equivalent position time constant of the closed velocity control loop.

The SINUMERIK 840D configuration for this machining centre does not include acceleration (torque) feed forward.

A velocity filter is included to damp the resonant frequencies in the closed position loop.

A velocity limitation in the form of saturation is included in the position loop.

Torque and current limitations in the form of saturation are included in the velocity loop.

The velocity controller does not include differential term (PI control).

Two additional filters are included in the velocity loop in order to get a filtering process with better time/frequency response. For example, Filter 1 can be configured as the PT₂ filter in the TNC 426 PB and filters 2–4 can be combined to get a band-rejection filter with better damping and frequency properties than the band-rejection filter in the TNC 426 PB.

The current controllers include integral term (see Motion controller section).

The sample time for the interpolator and position controller is 4 ms (cf. 3 ms).

The sample time for the velocity controller is 0.125 ms (cf. 0.6 ms).

The sample time for the current controller is 0.125 ms (cf. 0.2 ms).

Although the control algorithm is distributed in two different units (SINUMERIK 840D and the plug-in control units), the dynamics acting on the x- and y-axes are modelled as in the single-axis test rig case. The model for the rear bearing mounting has been updated to reproduce the fixed-supported bearing configuration.

The two-axis system is configured as shown in Figure 10(a). Linear and circular interpolation methods described in the Motion controller section are included in the interpolator to coordinate the movement of the axes.

Ford et al. ²⁵ derived algorithms for the real-time correction of time and spatial errors in a similar vertical machining centre configuration. Pre-calibrated geometric errors are included in the form of an error map that is used to correct the control movements over the working zone. Figure 10(b) represents the block diagram for the two-axis model. These errors are measured by using linear interferometers, straight edge, electronic levels, dial gauges, ballbars and granite squares.^25–27

The two-axis equations for the X–Y plane are derived from Postlethwaite ²⁶ considering only the 13 geometric components associated with x- and y-axes, thus

$E_x={\Lambda }_x(x)+{\Lambda }_x(y)+{\Lambda }_x(z)+{\phi }_y(x,y)D_z\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}+{\phi }_y(y,x)D_z+{\phi }_z(y,x)D_y+{\theta }_{\mathrm{xy}}(x,y)D_y$

(133)

$E_y={\Lambda }_y(x)+{\Lambda }_y(y)+{\Lambda }_y(z)+{\phi }_x(x,y)D_z\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}+{\phi }_x(y,x)D_z-{\phi }_z(x,y)D_x$

(134) The procedure presented by Ford et al. ²⁷ was used to isolate the load errors from the effects of the rigid body geometric errors. Table 3 contains the summary of the maximum values of the geometric and non-rigid errors measured from the actual machine.

Table 3.

Measurements of the geometric and non-rigid errors.

Axis	Component	Units	Geometric effect
X	Pitch: φ y(x)	µm/m	0.0484
	Yaw: φ z(x)	µm/m	(−5 to 2.8)*10−3
	Roll: φ x(x)	µm/m	(−7 to 6)*10−3
	Linear positioning: Λx(x)	µm	−0.05 to 9.15
	Horizontal straightness: Λy(x)	µm	−3.46 to 3.5
	Vertical straightness: Λz(x)	µm	−8.84 to −0.64
Y	Pitch: φ x(y)	µm/m	–0.0185
	Yaw: φ z(y)	µm/m	(−5 to 2.8)*10−3
	Roll: φ y(y)	µm/m	0 to 0.028
	Linear positioning: Λy(y)	µm	−0.14 to −13.92
	Horizontal straightness: Λx(y)	µm	−0.1 to −6.53
	Vertical straightness: Λz(y)	µm	−0.05 to −3.5
	XY squareness: θ(x,y)	µm/m	−26

Simulation approach for TLM models

Choice of simulation method

One of the aspects that lead to the study of the TLM method for the modelling of feed drives was the possibility of formulation of comprehensive models which would be implemented in real time. MATLAB^28,29 was used as the computational environment where high-level programming and visualisation functions were integrated for modelling, simulation and analysis of dynamic systems. The models were formulated in SIMULINK ³⁰ and then into a program or as block diagrams using the Graphical User Interface (GUI).

SIMULINK contains a large library of pre-defined blocks that supports the modelling of linear and non-linear systems in continuous time, sampled time, or a hybrid of the two. Systems can also be multi-rate, i.e. have different parts that are sampled or updated at different rates. SIMULINK features a tool called Real Time Workshop (RTW), which automatically generates C code from the SIMULINK models to produce platform-specific code.

Simulation method and development

Five SIMULINK multi-rate subsystem blocks, i.e. position controller, velocity controller, current controller, torsional loop and axial loop (the latter two including for mechanical vibrations and non-linearity’s such as backlash and friction forces) were implemented into the software for each axis. The variables interfacing the multi-rate subsystems were implemented in Data Stored Memory Blocks (DSM). A DSM defines a memory region for use by the data store read and data store write blocks. The feature gives access of the memory region to the different sub-systems in order to read from or write to at predetermined designed variable sample rates.

As an example to the approach, implementation of the single-axis model for the test rig in SIMULINK had to deal with the multi-rate requirements (see the Motion controller section): (a) position controller: v_ref at a rate t_p = 3 ms; (b) velocity controller: i_qref at a rate t_v = 0.6 ms and (c) current controller: e_dqref at a rate t_c = 0.2 ms.

In addition, each PWM signal is composed of seven e_dq voltages (switching states) calculated according to the expected currents to be induced in the motor. The duration of each e_dq voltage is specified in multiples of the propagation time for the torsional model (t_t). To accomplish this, t_t is made equal to the sampling period on the PWM signal (t_pwm). The propagation time on the axial model is a sub-multiple of the torsional propagation time as defined by the Castaneda method ⁹ which describes in detail the synchronisation of the torsional and axial models. These actions result in five multi-rate subsystems per axis which were implemented in software for multi-axis system application.

Model development

The modelling development was carried out according to the following methodology:

The model for the single axis test rig was built in SIMULINK in order to validate the modelling approach. The validation of the model was achieved by comparing simulated results with the experimental data recorded from the controller.

Following this step, the single and two axis models for the CNC machine were implemented in SIMULINK taking as a basis the TLM model for the test rig. The x- and y-axes models were validated against experimental data recorded from the controller and the axis drives. The controller's facilities included an integrated oscilloscope and the ability to input stimuli such as PRBS, jerk, and linear/circular interpolation. It also provided an additional capability to perform spectral analysis. The two axis model was validated against experimental data recorded in real time for a ball bar circular test and jerk limited responses (see section on Validation).

The models to include for the effect of geometric errors measured with specialised equipment such as laser interferometer, ball-bar, electronic precision levels, precision squares, etc. (Table 3). The simulated results were compared with the measured ones (see section on Validation).

Built-in functions on the CNC machine motion controller enabled spectral analysis to be carried out for the optimisation of the position, velocity and current closed loops for each axis. The control frequency characteristics were calculated by entering a PRBS at the set point of each control loop. Four frequency responses were measured, i.e. closed position controller loop, closed velocity controller loop, system and mechanical frequency response. The closed loop responses for the position and velocity controller were analysed and used for the optimisation of their parameters. The control system and mechanical frequency responses were used to set up the parameters for the filters that are used to damp the system resonant frequencies. The data measured were included in the system simulation and used also for the identification of axis drive modal parameters.

Other stimuli waveforms were entered to identify parameters suitable for the model such as step (velocity bandwidth overshooting check), triangular (backlash/ dead band), trapezoidal (position gain/coulomb friction) and random white noise/swept sine (resonant frequencies/damping coefficients).

The x-axis of the CNC machine was further modified in order to explore the possibility of a real-time implementation for the models. In this regard, the real-time workshop capability of MATLAB/SIMULINK was used to generate a real-time application targeting a RTI-1005 dSPACE environment. The RTI-1005 dSPACE platform was selected because it features the possibility for data logging, control and monitoring of systems in real time.

MATLAB/SIMULINK programs

SIMULINK programs compiled and developed

A bench mark MATLAB/SIMULINK program was compiled and applied to a previously simulated rotor shaft already in the public domain which was then used to evaluate the reduction of mathematical operations and mean square modelling errors when comparing TLM and MTLM techniques (Conclusions section).

Single Axis test rig model simulation (see section on TLM method applied to CNC machine tool feed drive)

Two axis CNC machine tool model simulation (see section on Implementation of two- axis TLM models to a CNC machine tool)

XY geometric calculation in a SIMULINK block. Developed for insertion of measured geometric errors at the machine into the simulated two axis model (see sections on Implementation of two- axis TLM models to a CNC machine tool and Validation).

MATLAB programs compiled and developed for the following

Structure of data for the pulse propagation on the axial model (see Screw shaft axial model section)

Parameters and initialisation code for the torsional loop subsystem (see Screw shaft torsional model section)

The PWM generating function (see sections on Inverter and motor (electrical) and Motor mechanical equations and coupling)

Calculation and nut monitoring blocks (see Screw shaft, nut and table section)

Parameters and initialisation code for the axial loop subsystem (see Screw shaft axial model section)

Identification of drive modal parameters using the Continuous Wavelet Transform (CWT) for improved performance when compared to other available analysis techniques

Automatic modal parameter identification using the CWT

Reference signal profiles swept sign, sinusoidal, white noise, jerk, linear interpolation and circular interpolation (for application comparison)

Axis block parameters and initialisation code (see section on Motion controller)

Calculation of polynomial coefficients for the geometric errors (see section on Implementation of two- axis TLM models to a CNC machine tool)

Calculation of the dilation parameters associated with CWT

Identification of damping factors and resonant frequencies to be used for axis control system frequency response.

The complexity and detail required for this section is beyond the scope of this publication. To address this issue, we have put all the detail for drive data, measurement methodology, parameter determination, coefficient analysis and a full description of relevant MATLAB/SIMULINK programs online. Most of the detail being contained within the appendices of that report. ⁹

Validation

Step velocity stimuli were tested for 500 mm/min (see Figure 11a) and 1000 mm/min where the simulated velocity and current responses (red line) match the machine measured values (green line). A 3% error was achieved in the transient area (0–0.02 s). The % error increases to 5% when the motor reaches constant velocity (Figure 11b).

The jerk-limited response was used to verify the performance of the axis drive during acceleration, deceleration and operation at constant velocity. For this purpose, position and velocity responses were measured for different displacements (10, 20, 100, 200, 400, 400 mm) and corresponding feed rates (500, 1000, 5000, 10,000, 20,000 and 40,000 mm/min). Figure 11(c) shows the comparison between measured and simulated position following error for a displacement of 200 mm and a feed rate of 10,000 mm/min. The position versus time measured at the machine and simulated responses matched very well with no significant error but the following error differed minimally by 2 µm at its peak in the acceleration zone and 8 µm at its peak in the deceleration zone. Figure 11(d) shows the comparison between the measured (green line) and simulated (red line) actual velocity signals over the zoom acceleration zone.

The same reference position and velocity feed-forward signals (radius 150 mm, feed rate of 1000 mm/min) were used for the circular interpolation in the X-Y plane at the machine and for the validation of the x- and y- axis TLM model responses (i.e. response to sinusoidal and cosinusoidal demands).

The two-axis model was simulated and the model position responses were introduced to produce the trace of x against y. The difference between the trace and the circle is illustrated in Figure 11(e) in the form of a circular plot analogous to that measured using a double ball bar.

The same sinusoidal stimuli used for the simulation was applied at the machine and the position outputs measured using a Renishaw ball bar system (Figure 11f). An angular overshoot of 180° before and after data capture for a two cycle 360° operation was utilised to minimise uncertainty of test synchronisation.

Figure 11(e) plot has an oval shape as a result of the squareness error (−26 µm/m) which is consistent with Figure 11(f). The TLM model shows that the introduction of the geometric errors generates a progressive error deviation when the machine worktable prescribes a circle. The simulated results do not match the machine measurements on the arcs described in the (45–90°) and (135–270°) zones but with sensitivity analysis this could be further improved. The difference may be due to divergences on the straightness measurements for the x-axis. The reversal spikes match closely on both the x- and y-axes.

Figure 11.

Validation responses comparing TLM model with measured machine. (a) Step velocity response, (b) reference current response to step, (c) position error response jerk-limited input, (d) velocity response to jerk limited input, (e) simulated ball bar trace and (f) measured machine ball bar trace.

Conclusions

A new MTLM stub was developed which improves the convergence and computational processing speed of the original stub algorithm. A reduction of 40% in the number of mathematical operations and almost 35% reduction in the mean square modelling error figures were achieved over the original TLM stub model.

Comprehensive transmission line models for the elements of a typical arrangement of a CNC feed drive have been described. All non-linear functions including geometric and load errors have been included as calculated and identified by the measurements undertaken. Equipment such as laser interferometer, ball bar, electronic levels, signal acquisition systems and others were used to obtain parameter data.

A TLM model for a CNC single- axis feed drive including a digital controller has been developed and then extended to a two-axis drive of a machining centre including geometric and load errors, thus forming the basis for universal CNC machine tool models.

The simulation of the single-axis and two-axis models to various feed rates and displacements including linear and circular interpolation matches well in comparison with the measured response of the machines under study (see Validation section).

The application of the MTLM transform and the torsional/axial model’s synchronisation approach to the modelling of single- and two-axis drives led to a real time implementation of the feed drive models. Subsequent work not reported in this paper has shown the feasibility of developing the model for implementation of the two- and three-axis models running in real time.

The x-axis of the machine was modified later in order to explore the possibility of a real-time implementation for the models. In this regard, the real time workshop capability of MATLAB/ SIMULINK was used to generate a real-time application targeting a RTI-1005 dSPACE environment. The RTI-1005 dSPACE platform was selected because it features the possibility for data logging, control and monitoring of systems in real time.