Journal of Dynamic Systems, Measurement, and Control

Moshe Shoham

Journal of Dynamic Systems, Measurement, and Control

2000, Journal of Dynamic Systems, Measurement, and Control

Abstract

Transactions of the ASME® Dynamic Systems and Control Division Technical Editor, A. GALIP ULSOY Past Editors, W. BOOK M. TOMIZUKA DM AUSLANDER JL SHEARER KN REID Y. TAKAHASHI MJ RABINS Associate Technical Editors, Y. CHAIT „2002… F. CONRAD „2001… E. FAHRENTHOLD „2001… S. FASSOIS „2001… Y. HURMUZLU „2001… R. LANGARI „2002… E. MISAWA „2002… S. NAIR „2001… N. OLGAC „2001… C. RAHN „2002… S. SIVASHANKAR „2002… J. TU „2002… P. VOULGARIS „2001… Book Review Editor A. GALIP ULSOY Executive ...

Figures (383)

Fig. 1 The rotational inverted pendulum system

Now, it is obvious that the right-hand side of the upper equation in system (34) does not depend on 6,. Following the approach of case 3, if the condition such that the right-hand side of the upper block in the motion equations would not depend on the time derivative of the state variable of the bottom block. Since x=3}—(C,/K,)6,, substitut- ing j' from system (31), yields

Table 1 Parameters of the rotational inverted pendulum system

Fig. 2 Hardware setup configuration of the pendulum system

Further, we will focus on the performance of the inverted pen- dulum system using our previously developed sliding mode con- trollers. Two case studies of the control objectives: case 1, stabi- lizing the pendulum at 6,=0 with @)=0, and case 2, stabilizing both the pendulum and the rotating base with respect to the equi- librium point 0)= 0,=0, will be presented.

Fig. 4 Simulation results by SMC for case 1: stabilizing the pendulum and it is fixed for other experimental results in the later figures. The control law utilizing a continuous approximation by a sinu- soidal function is designed as

Fig. 6 Experimental results by SMC: sloshing water (case 1) Fig. 5 Experimental results by SMC: no weight (case 1)

Fig. 7 Experimental results by SMC: metal bolts (case 1)

Fig. 8 Simulation results by SMC for case 2: stabilizing both the pendulum and the base

Fig. 9 Experimental results by SMC with small 6: no weight (case 2)

Fig. 10 Experimental results by SMC with larger 6: no weight (case 2)

Fig. 11 Experimental results by SMC: sloshing water (case 2)

Fig. 12 Experimental results by SMC: metal bolts (case 2)

Application of the Passivity Theorem. Using the notion of multipliers [2], the feedback systems shown in Fig. | are equiva- lent from the point of view of stability. The passivity theorem states that if G is passive and H is strictly passive then r eL,=yeL,. Using Lemma 1 and the passivity theorem, the original system is stable if either H= Y, TH or H= HY, is strictly passive. A possible design strategy would select H,,, a strictly Therefore G(s) is passive if and only if w<* When «= p*, the summation over the vibration modes vanishes (all vibration modes become unobservable) and the output p,,+=C,i7, becomes propor- tional to the rigid body modal rate. It 1s instructive to consider a small negative rate feedback u(t) = — ep,,. Using a simple eigen- value perturbation argument, it is readily shown that the closed- loop eigenvalues of each vibration mode are given by —€Y Cr lratja- All modes are stabilized when u<p* (Y, >0) and all of them are destabilized when p> (Y ,<0) as- suming controllability and observability, i.e., ¢.,#0.

Fig. 3 Torsional control system apparatus (reprinted courtesy of Educational Control Products)

The noncollocated degree of freedom is taken to be p,.= 03= 9 +[1 1]q, so that C,=1, C,=[11]. The collocated degree o' freedom is p.o= 9, and p,=wA,;+(1—) 4. If g-=6,, de) = 92.— 9, and g,.= 03— 9» are selected as gen- eralized coordinates and f(t) as the motor torque, the matrices defined in Eq. (2) are given by

Fig. 5 Parameter estimates (4*=0.973, 4=0.96)

Fig. 6 Tracking errors and parameter estimates (u*=0.567, p=0.55) Table 2 Measured values of 4, at stability-boundary

pole locations p, are defined by p.= Bp, (see Fig. 1), or if dis- crete time control is used p.=e P47, where AT is the sampling time. This way of scaling maintains the damping of the poles, while > changing the speed of the response.

Fig. 2 Scaling design for PD control of a double integrator plant From Fig. 2 it is clear that the performance can be improved by increasing 6 outside the interval (0,8);,) where the control does

Theorem 4.1 is illustrated schematically for a second order ex- ample in Fig. 3. Theorem 4.1 does not say anything about how to obtain the matrix P. Any choice of P that satisfies Eq. (8) for a =1 can be used; it then follows that there exists an ayin<1.

Table 1 Comparison of design methods, full state feedback of a number of numerical optimization techniques can be used to minimize S(P,) for a fixed 8, and thereby find a P such that

Fig. 4 Simulation results for x(0)=[0 0.71 0.71]? and B=0.79, 2.29, and 3.23 conditions is a limited number of points in state space, it is pos- sible to obtain 6,,, quite accurately. The accuracy is determined by the size of the change in the parameter 8 used when searching for Boot -

with plant state x e %", integral state x; = 9%, and controlled output ze. The AW is shown in block diagram form in Fig. 5, where q denotes the unit delay operator. This structure is very typical for many AW methods; what is special about this method is the spe- cial choice of the AW gain that results in the integral state becom- ing a function of the state x during saturation.

Table 2 Comparison of design methods, full state feedback plus integral control

Fig. 10 Experimental and simulation result for full state feed- back control 9 Conclusions Finally, the stability based method is applied to the case of full state feedback plus integral control for the same plant as in the previous example. In this case it turns out that the limiting factor is not the saturation, but the unmodeled dynamics. This means that the design parameter £ has to be limited to avoid exiting the

Fig. 2 Schematic of a multilimbed robot Fig. 1 Simplified Cartesian computed torque (SSCT) control

ee ee ee ee eee eee Oe eee ee The SCCT control scheme is applied to a multilimbed robot, such as the robot shown in Fig. 2 as follows. The various limbs of the robot are viewed as manipulators, each doing a different task. Some limbs support and position the robot body, some carry and manipulate objects, and some apply forces to the environment. Each limb is under SCCT control and the desired reference inputs are appropriately selected for its task. For example, for free limbs that move in space, the controller shown in Figure 1 can be ap- plied with X,,, and X,,, being the desired and actual position/ orientation of the tip of the limb with respect to the body. For limbs that are in contact with the ground and are used to position the body, the variant of SCCT control shown in Fig. 3 can be used. In this case, Xq., and X,,; are the body’s desired and actual

Table 1 Model parameters of LIBRA’s joints

Fig. 3 Block diagram for the control of a multilimbed robot

5 Experimental Results of SCCT Control Applied to One Limb

Fig. 6 Comparison of experimental system and simplified model Experiments were performed to compare the Jacobian transpose and SCCT control schemes for simple endpoint under position control of one limb. Figure 7 compares the two controllers for a step input. The tip of the limb is commanded to move 2 cm in the Y direction: from point A to point B in the figure. The SCCT controller is able to hold the limb’s X position essentially constant

Fig. 7 Step response of two Cartesian controllers

Fig. 8 Experimental path tracking of the two Cartesian con- trollers

Fig. 9 Stiffness selection in Cartesian space

Fig. 10 One complete step in a LIBRA climbing gait

Note that, in the basic backstepping procedure, at Step n one determines the actual control as aS Eee ee ee eee lee ee To apply the second order variable structure control strategy within the backstepping framework, the sliding manifold

The two dynamical Eqs. (12) complete the transformation from the state [x1(1),....x,(D ]7 to the state [21 (1),.--5Z,—1(t),e1(t),e2(t)]’.. With implicit symbol defini- tions, system (12) can be rewritten as The idea is now that of steering both e,(ft) and the unmeasurable (ft) to zero (to attain the sliding regime on S(t)=S(t)=0) ina finite time, that is, a second order sliding mode control problem has to be solved. In Bartolini et al. [26], it has been proved that, by analogy to the well-known solution to the time optimal control problem, the control u(t) can be chosen as a bang-bang control switching be- tween two values — Uy, and +Uy,,. The classical switching logic for a double integrator (H(e,(t), e2(t))=0, Bj, =B,=1) is

J, J ace the load and motor inertias, F,,F,, are the load and motor damping constants, k is the joint stiffness, and u is the input torque applied to the motor shaft. Assume that J,,=J,,(x) as it happens in practical situation in which the motor is driven by the input voltage through a nonlinear element the characteristic of which is a function of the system state. System (25) can be rewrit- ten as

VxeB,, V@eBj, where B, and Bg are neighborhood of the equilibrium point in the x-space and 6-space, respectively. B(x) has been replaced by B(x,v), that is a linear function in the v variable (it represents the dependence from the @ parameters in Eq. (2), here is substituted with its estimate as it appears as coef- ficient in the expression of u).

Fig. 2 Output trajectory in the proposed procedure

Fig. 5 S trajectory in the proposed procedure

Fig. 4 u trajectory in the proposed procedure

Fig. 6 0 trajectories in the proposed procedure

where q, g, and g are, respectively, the n <1 joint position, ve- locity, and acceleration vectors, u represents the n X 1 torque vec- tor generated by actuators, M(q) is the symmetric positive defi- nite generalized inertia matrix, C(q,q) is the force (torque) vector resulting from Coriolis and centripetal accelerations, G(q) is the generalized gravitational force vector, and d(q,q) denotes the dis- turbance. Define a 2n-dimensional state vector x as

This completes the proof of the following theorem: ee - Theorem: For a robot model Eq. (1) subject to sliding mode based repetitive learning control, which is accounted for by Eq. (19), the sliding surface S and the tracking error e are uniformly bounded if both gain matrices K and Q in the reaching law are adequately chosen. Furthermore, having learned a number of cycles, the ultimate tracking error e is bounded by [t is always feasible to adequately choose K, A, T’, and oy such that R is positive definite. Therefore, one can prescribe values of positive constants o,, 7,, and og to satisfy

Fig. 1 Schematic diagram of direct-drive robot As shown in Fig. 1, this study constructs a three-axis R— 6 —Z direct-drive robot manipulator, where the first link is driven by a NSK Megatorque motor [14], the second by a NSK Mega- thrust motor, and the third by an electrothrust motor together with ball screw.

The speed and acceleration profiles of the end-effector are pre- scribed as shown in Figs. 4(a) and 4(b), respectively. Figure 5 depicts a spatial circular trajectory to be tracked. In addition, a planar square trajectory on the X—Z plane, will also be carried out in experiments. Figure 6 depicts the control block diagram. To implement the present method, Eqs. (11), (13), and (14) are re- written to become the discretized form:

Fig. 6 Control block diagram of the present method

Fig. 4 (a) Speed profile of end-effector and (b) corresponding acceleration profile

Fig. 8 Position errors in Z-axis using (a) sliding mode control and (b) the present method Fig. 9 Sliding variables in R-axis using (a) sliding mode con- trol and (b) the present method (1) Spatial Circular Trajectory Tracking. In contrast to Fig. 8(a), tracking errors in Fig. 8(b) progressively reduce to a smaller

Fig. 13 Tracking result of the first time learning in X-Z plane Fig. 10 Tracking result of the first time learning in Y-Z plane

Fig. 12 Position errors in Z-axis using (a) sliding mode con- trol and (b) the present method Fig. 11 Tracking result of the fifth time learning in Y-Z plane

Fig. 14 Tracking result of the fifth time learning in X-Z plane

and where use was made of the definition of E(x) and of the fact that po;#0, ie J,. Next consider the matrices [[9(x) Ao(x)] and [T,(x) A,(x)] defined as

We are now able to establish the following theorem. Eee eS ee ee eee ee Theorem 4.1. The DRLD problem via the restricted state feed- back law (2.3) is solvable if and only if there exists a set of nonnegative integers {Y1,s,-.-,Yp,s+ Such that the following condi- tions simultaneously hold:

Example 5.1. Consider the nonsquare nonlinear analytic system of the form (2.1) with whose input-output description consists of the following set of differential equations: In order to demonstrate the efficiency of our approach, the fol- lowing illustrative examples are presented:

Example 5.2. Consider the nonsquare nonlinear analytic system of the form (2.1) with Finally, as it can be easily checked, for the choices (y\,;,¥2,) =(0,1) or (1,0) or (1,1) or (0,2) or (2,0), Theorem 4.1 cannot be satisfied.

Hence, the appropriate restricted state feedback is characterized by where fo; .i€ J, are arbitrary nonzero real parameters. It is easy to see that, application of the above restricted state feedback to the system at hand yields an I/O linearized and simultaneously decou- pled and disturbance localized system, whose input-output de- scription has the form

Therefore, the appropriate restricted state feedback is given by where fo;, 1¢Jy are arbitrary nonzero real parameters. The closed-loop system has the input-output description It can be easily checked that for the remaining choices of the pair (Vis+¥2s), Theorem 4.1 cannot be satisfied.

where f9;,1€J2 are arbitrary nonzero real parameters. The closed-loop system has the input-output description Therefore, the appropriate restricted state feedback is given by

Space station dynamics: Journal of Dynamic Systems, Measurement, and Conirol NE Rg OE Ley © Example 5.4. The nonlinear equations of motion of the space station rotating in circular orbit with constant orbital rate, in terms of components along the body-fixed control axes can be written as [28]

Assuming a constant orbital rate and differentiating (5.5), we obtain Solving (5.3) with respect to the vector [w, w w3]|, we obtait

where, g;,92.q3 are the body-axis components of the CMG mo- mentum. Pa a Sy aa ee Oe ee Ee 2) 317T .— .44¢.°. Attitude kinematics (2-3-1 body-axis sequence):

Therefore, Theorem 4.1 is satisfied for (y1,5,¥2,;)=(0,0) and the problem is solvable. Moreover, we have

Fig. 1 Sensitivity curves relative to normalized frequency for ZV, ZVD, and El input shapers (@mogei=1 rad/sec, Cmoget= 0-01, ty=4a sec).

Fig. 2 Shaper lengths for optimal shaper designs assuming uniform (top) and Gaussian (bottom) variation in the natural frequency. The ZVD and El shaper lengths (which are equal) are also shown for comparison.

Fig. 3 Comparison of the 5 percent insensitivity W; levels as- suming for optimal shaper designs uniform (top) and Gaussian (bottom) variation in the natural frequency.

Fig. 4 Comparison of the expected residual vibration (J) assuming for optimal shaper designs uniform (left) and Gaussian (right) variation in the natural frequency. For smaller parameter uncertainties, the ZVD is more optimal; and as the parameter uncertainty increases, the El shaper becomes closer to optimal.

Fig. 5 Comparison of the magnitude of shaper zeros assuming for optimal shaper designs uniform (left) and Gaussian (right) variation in the natural frequency.

Fig. 6 Shaper length of minimal expected residual vibration shapers assuming uniform and Gaussian variation in both fre- quency and damping. For lightly damped systems, optimal shapers yield shorter shaper lengths than those obtained from traditional ZVD and El shaper designs.

Fig. 7 Expected residual vibration (J) assuming for optimal shaper designs uniform (top) and Gaussian (bottom) variation in both frequency and damping. Optimal shapers yield lower expected residual vibration than those of ZVD and EI designs.

Fig. 9 Shaper lengths for different uniform frequency uncer- tainties. The optimal shapers yield shorter shaper lengths com- pared with traditional ZV shaper designs. As in the comparison of three-impulse optimal shapers and tra- ditional ZVD and EI shapers, we see that the optimal shapers lead to faster maneuvers while yielding lower levels of expected re- sidual vibration, as shown in Figs. 9 and 10. A uniform frequency

Fig. 8 Angle deviation of optimal shaper zeros from system poles for uniform (top) and Gaussian (bottom) variation in fre- quency and damping. The angle deviation is negative for the optimal shaper zeros meaning that their location is shifted to- ward a line of lower damping.

Fig. 10 Expected residual vibration for different uniform un- certainty ranges in frequency. The optimal shapers yield lower expected residual vibration compared with the traditional ZV shaper designs. uncertainty size of 0.1 means w is uniformly distributed from 0.9@ moder tO 1.1@mode, and similarly for other frequency uncer- tainty sizes.

Fig. 2 Material deposition from mass source to the molten globule

Fig. 3 Numerical simulation of thermal material transfer

Fig. 4 Molten globule modeling through a spherical lens primitive However, when solidification occurs under nonequilibrium con- ditions, i.e., during surface oscillation, or under the influence of the mass and/or heat source and a shielding gas stream, the geom- etry of the elementary deposition field deviates from the simple shape of Fig. 4. The same is true when material is deposited on a nonflat substrate (i.e., adjacent to previously laid beads) or not horizontal to gravity. Figure 5(a) shows such a unit deposition distribution identified in gas metal arc welding (GMAW) with a consumable stainless steel electrode over a similar substrate. Fig- ure 5(b) illustrates how such a deposition field Z can be composed by a number n of elementary spherical lens primitives Z;, each defined by parameters u;,U;,W;.1;:

Fig. 5 Material deposition field Z in GMAW: (a) composite deposition distribution Z; (b) decomposition into spherical primitives > Z;

Fig. 6 Distributed-parameter model of surface geometry z by superposition of deposition fields Z Equation (9) updates the time-dependent, arbitrary-shaped deposition distribution Z, by assessing the deviation field «=z —z between the actual measurements z and the predicted geom- etry z. This distributed-parameter identification can be performed in real time, by periodic adjustment of Z every Dt, with only a limited number of local geometry data at locations (y,y) needed to

Equation (8) actually represents a locally linearized, time-varying (LTV) description of the nonlinear phenomena composing mate- rial deposition, as already discussed. The nonstationary lineariza- tion and the respective parameter uncertainty of the mass transfer nonlinearities are reflected in the time variation of the deposition distribution Z. Thus, rather than anticipating an exacting model of these complex effects on physical grounds, an approximation of the composite deposition field Z Eq. (5) will be updated through experimental identification of its parameters (u;,v;,w;.r;). If laboratory data on a part of the surface geometry z(x,y;f) are obtained over a time period Dt, the deposition field Z can be identified by incremental deconvolution of Eq. (8), based on the applied material feed F(X,Y;7) from the mass source along its motion:

In principle, given geometric measurements z at 4n locations (y,y) on the part surface, and the respective predictions z of the deposition model from Eq. (8), the resulting deviation values can be used in Eq. (10) to form a linear system of 4n equations, to be solved for the deposition parameter changes (Du; ,Dv;,Dw;,Dr;),i=1...n. Notice that in order to obtain such a well-conditioned system of equations for robust parameter identification, sufficiently rich data must be obtained at respective 4n positions (y,), located close and scattered around the current material deposition region. In the vicinity of the mass source (X,Y), however, real-time measurements, e.g., by mechanical pro- filometry or optical scanner sensing, are technically problematic because of the material deposition interference with the sensor. Moreover, the precision of the geometry model is challenged by the variable source influences, due to the dominating surface ten- sion assumption. For this reason, in-process geometry measure- ments at m locations are conducted on the solidified surface at a safe distance behind the mass source, as explained below. On the basis of these m measurements, the 4n deposition parameter changes are evaluated by a least-squares identification technique (Astrom and Wittenmark [24]). This method minimizes a qua- dratic index *e(t) of the deposition error ¢ values at the m mea- surement locations, analogous to that of Eq. (6). The resulting corrections Du;,Dv;,Dw;,Dr;,i=1...n, are subsequently used every time period Dt to update the deposition field param- eters u;,U;,W;,r;, to their new values u;+Du;,v;+Dv;,w; +Dw;,,r;+Dr;. identify the parameters (u; ,v; ,w;,7;) of the composite deposition field approximation Z, Eq. (5). For this purpose, Eq. (9) can be rewritten over the full process duration ¢, as a linearized expres- sion of the deviation field ¢ sensitivity on the parameter variations Dy. Dn. Dw. Dr. j=1.. 2 -n-

Fig. 7 Experimental station for GMAW thermal material depo- sition with laser stripe optical scanner The geometry of the deposited material profile is measured in- process, along a line segment across the direction of motion and trailing behind the mass source (Fig. 8). This is obtained by a 3D optical scanning system based on a sweeping laser stripe. It con- sists of a 20 mW laser diode, which through a cylindrical lens projects a bright light stripe on the prototype surface. The laser stripe is shed across the deposited bead and it follows the GMAW torch at a distance I=25.4mm, separated by an incandescent cloth partition to minimize light interference from the arc. The deflections of this line as it sweeps the deposited material reflect the geometric profile of the latter, and are monitored by an optical CCD camera with band-pass filtering at the laser wavelength (680 nm). The optical image is monitored in real time by the system computer, through a frame grabber and image processing soft- ware, and is also recorded in standard video format for off-line analysis. This software includes algorithms for the CCD camera elevation angle correction and gray-scale interpolation of the height z at the m measurement locations O; ,i=1...m along the laser stripe, with an accuracy of 0.1 mm. Note that the distance I

Fig.8 Arrangement of laser stripe scanner for geometric mea- surement of bead profiles

Fig. 9 Open-loop bead width and height responses to torch velocity changes: zig-zag line: experimental measurements; long-dashed line: analytic model predictions

Fig. 10 Block diagram of distributed parameter geometry controller for material deposition

Fig.1 Stribeck curve. The duty parameter is defined as p;v/F, where yp; is the lubricant viscosity, v is the sliding velocity, anc F,, is the normal force per unit contact area.

Fig. 2 A simplified condition of polymer-to-metal contact with lubricant. m, is the average separation between the polymer and metal surfaces with maximum My and minimum M,,, A, is the local deformation on the polymer surface with maximum Ay and minimum 4A,, .

Fig. 3 Schematic of a hydraulic actuator, the seals are all made of polymers and both sides of the piston are filled with hydraulic fluid. Lubricated polymer friction exists at all seal-to- metal contacts.

Fig. 4 Comparison of the deterministic model (dashed line) and the measured data (solid line)

Fig. 5 Comparison of the overall stochastic model (dashed line) and the measured data (solid line)

Fig. 6 The estimated normalized autocorrelation function of the prediction error of the stochastic model. The two parallel dashed lines denote 95% confidence interval.

Fig. 8 The estimated normalized autocorrelation function of the prediction error of the extended Kalman filter after 1 hour and 32 minutes of reciprocating sliding. The two parallel dashed lines denote 95% confidence interval. Fig. 7 A close-up of the extended Kalman filter’s estimates (marked by *) and the measured data (solid line) after 1 hour and 32 minutes of reciprocating sliding

Fig. 9 A close-up of the extended Kalman filter’s estimates (marked by *) and the measured data (solid line) after 8 hours of reciprocating sliding

5 Conclusion Fig. 12 The trajectory of parameter a, after 8 hours of recip- rocating sliding

Fig. 11 The trajectory of parameter a, after 1 hour and 32 min- utes of reciprocating sliding Fig. 10 The estimated normalized autocorrelation function of the prediction error of the extended Kalman filter after 8 hours of reciprocating sliding. The two parallel dashed lines denote 95% confidence interval.

Fig. 13 The trajectory of parameter G3 after 1 hour and 32 minutes of reciprocating sliding

Fig. 14 The trajectory of parameter G3 after 8 hours of recip- rocating sliding

Fig. 2 Flow chart showing the battery of tests for identifica- tion and characterization of underlying process dynamics Fig.1 Poincare section plots of force and vibration signals—collected at cutting speed V=130 ft/min, feed f=0.0088 in./rev—capturing the variations in machining dynamics with flank wear. Fresh, partially worn and fully worn refer, respectively, to h,=0.000, 0.0060, and 0.0175 in.

Fig. 3 Methodology of fractal estimation consisting of an offline training phase and an online estimation phase 3.1 Experimentation. The experiments were conducted on a 20 HP LeBlond heavy duty lathe. The workpieces were made of 36in.x é7in. SAE 6150 Cr-V steel, and the tool inserts were uncoated carbide grade K68 with geometric specification SPG- 422. Three different online sensors were used—(i) a three-axis Kistler Z3392/b piezoelectric dynamometer, placed underneath

Fig. 4 Design of experiments: Solid and dashed outlines of circle indicate whether that the exemplar patterns corresponding to that design point are used for training or testing

Fig. 5 Representative plots showing (a) the variation of aver- age fractal dimension with d; for two experimental runs con- ducted using a fresh tool: V=160 ft/min, f=0.0136 in/rev. (b) the effect of signal separation on the fractal dimension esti- mates: — corresponds to the separated signal, and —--—-—-— corresponds to the nonseparated signal. As a result of signal separation, the change of slope of the graph becomes more pronounced. This reduces the uncertainty in deciding the linear portion of the graph thus rendering the computation of the slope of the linear portion of the graph, and hence the fractal dimension estimates more accurate.

Fig. 6 Architecture of the recurrent neural network Table 1 Sample input patterns Table 2. Neural network training details

Fig. 7 Influence of time of cutting, and hence the extent of flank wear, on the estimation error

Fig. 8 (a) Comparison of interpolates of the flank wear esti- mates (— — —) and the actual measured flank wear (—). The experiment was conducted at V=160 ft/min and f=0.0136 in/rev. (b) An empirical q—q plot showing the equivalence of the distributions of the estimation errors corresponding to the exemplar patterns used for training the neural network, and those corresponding to the testing patterns.

Fig. 1 Flowchart for the FRF-based identification scheme The overall procedure for joint parameter identification is sche- matically depicted in Fig. 1. The identification of joint parameters requires an FEM model of the structure without joints and experi- mental measurement of FRF’s at a few selected points. FRF’s at necessary points, particularly at joint locations, are estimated by using an indirect estimation technique described in Sec. 2.1. The joint parameter identification procedure is described in Sec. 2.2.

and W;(w) is a weighting function to take into account the reli- ability of identification with respect to frequency. Assuming, without loss of generality, the bearings have con- stant stiffness and damping coefficients over a certain range of frequency, it will be useful to make use of the following weighted average for the bearing parameters, i.e., Since the above identification formula has a linear relationship with the indirectly estimated FRF’s, it may be useful to use //(w) as the weighting function for Eqs. (11) and (12) or for selecting the frequency range appropriate for identification. In addition, magnitudes of the FRF’s involved in the matrices can reveal the relative importance of the results as a function of frequency. Therefore, the magnitudes of some FRF’s related to the identifi-

Fig. 2 Configuration of the high speed spindle with cutter and toolholder 3.1.1 Experimental Setup and Measurement of FRF’s. The spindle system consists of a 6 in. cutter connected to the spindle by a size 50 taper toolholder. The spindle is a horizontally mounted high speed spindle with the maximum rotational speed of 10,000 rpm and is driven by a 100 HP DC motor through a 3 in. wide belt. The spindle is equipped with two sets of tapered

Fig. 3 Experiment 1: (a) index functions; (b) measured FRF at the rear end; (c) measured FRF at the front 3.1.2 Identification of Stiffness Coefficients. The selection of the frequency range for identification is a two step process. The first selection is based on a broader range of frequency and is used in the identification program. This initial range is determined by cluding the majority of the main lobes of the frequency re- sponse. For experiment 1, this range was chosen to be from 650 to 300 Hz based on Fig. 3 (b-c). The identification produces an identified stiffness value at each frequency in the selected fre- quency range. The second frequency range is determined by a

Table 2 Identified stiffness of the high speed mill spindle combination of the index functions and the magnitude of the FRE. Once the second frequency range is determined, a weighted aver- age based upon the index functions provides a single, averaged stiffness value for each bearing to be identified. The second quency range selected for experiment 1 was 760 to 790 Hz d fre- eter- mined from Fig. 3 (a—c). To verify that the measured FRF’s are acceptable, coherence functions were obtained. The coherence functions for experiments 1 and 4 are shown in Fig. 4. Since experiments 5 through 7 were conducted with the changed p hysi- cal system, the first frequency range selected for the identification was changed accordingly to include the majority of the main lobe measured for each experiment. The second frequency range se- lected corresponded to the higher index values associated with the main lobe.

Fig. 5 Comparison of measured FRF’s and estimated FRF’s with transverse damping for experiment 2: (a) rear; (b) front The stiffness values for all the experiments are provided in Table 2. Since the front and rear bearings used in this milling machine are two different types, the stiffness parameters are not

Fig. 4 Coherence functions in FRF measurement: (a) experi- ment 1, front; (b) experiment 1, rear; (c) experiment 4, front; (d) experiment 4, rear

Fig. 7 Comparison of measured FRF’s and estimated FRF’s with transverse and tilting damping for experiment 2: (a) rear; (b) front Fig. 6 First three mode shapes of the horizontal mill spindle

3.1.4 Tilting Damping Identification. Observation of experi- mental results shows that the third mode is overdamped. How- ever, since the bearings (located at 0.2083 and 0.5458 m) are near the nodal points for the third mode as shown in Fig. 6, transverse damping alone will not be able to correctly model the damping

Table 3 Identified damping of the high speed mill spindle

Fig. 8 Comparison of measured FRF’s and estimated FRF’s with transverse and tilting damping for experiment 1: (a) rear; (b) front

Fig. 10 MAZAK spindle experiment: (a) index functions; (b) measured FRF on the toolholder directly below the spindle; (c) measured FRF at the tip of the extension tool

Table 4 Theoretical and identified stiffness for the MAZAK spindle

Table 5 Natural frequency prediction and validation for the MAZAK spindle Nomenclature

In order to rearrange the system matrices into measured and un- measured DOFs, an identity relation is introduced into Eq. (A2)

Fig. 1 Probe type sensor and CCS: (a) using four probe sen- sor; (b) using the CCS

Fig. 4 Capacitance at the eccentric rotor position Wiodeling OF the KouNndness Krrors Of the Kotor and the Sensor. The roundness error of the rotor deteriorates the distor- tion of the line, as shown in Fig. 4. If the roundness error is h (function of 6), the small variation Au from the nominal straight u,, can be expressed as Fig. 3 Capacitance at the concentric rotor position

Fig.5 The effect of the sensor size ¢ on the error amplification factor e,,/h,, with various odd harmonic numbers of m

Fig. 7 Test rotors with harmonic roundness errors: (a) rotor 1; (b) rotor 2

Fig. 6 Measurement errors of two probe sensors, four probe sensors, and the CCS (dia. of rotor 50 mm gap between sensor and rotor 0.5 mm magnitude of roundness error 0.001 mm ra- dius of revolution 0.003 mm)

Fig. 8 Magnitudes of harmonic errors of each testing rotor

Fig. 10 Simulation of measuring orbits of test rotors with vari- ous sensor types (dia. of test rotor 50 mm depth of harmonic error 0.05 mm, radii of revolution 0.02 mm, and 0.03 mm). (a) Rotor 1 with third harmonic error; (b) rotor 2 with fourth har- monic error. Two test rotors are shown in Fig. 7 and the experimental setup is built as shown in Fig. 9(a). The experimental setup consists of speed controlled DC motor, hydrostatic bearing, oil-less bearing, spindle with test rotor and sensor housing with four probe sensors, and the CCS. The probe sensor used is a capacitive sensor, PX405HA manufactured by LION PRECISION. The photo of the

Fig.9 Experimental setup: (a) schematics; (b) sensor housing

Fig. 11 Experiments of measuring orbits of test rotors with various sensor types (dia. of test rotor 50 mm, depth of har- monic error 0.05 mm, radii of revolution 0.02 mm and 0.03 mm). (a) Rotor 1 with third harmonic error; (b) rotor 2 with fourth harmonic error.

Fig. 13 Mounting errors of probe sensor small rotation of the sensor coordinates from the reference coor- dinates and by not heading toward the reference origin. As the latter exists at =O, the effect of this error is expressed as small distortion of the sensor.

Fig. 12 Mounting errors of the CCS: (a) misalignment; (b) angle error

Fig. 1 The unperturbed system in the roll manifold

Fig.5 Roll dynamics with parameter variations under different controllers: nonlinear feedback (solid), linear feedback (dash— dot), and open loop (dashed) Without parameter variations, the peak control effort using non- linear feedback control is uy,,,=1.45 for IC1 and u,,,,=2.24 for IC2, and for the linear feedback control that value is uy,,,=3.40 for IC2. With parameter variations under nonlinear control, upax =3.88 for IC1 and w,,,,=8.05 for IC2. These peak control efforts usually occur during the initial transient period and settle down to smaller levels soon after. The peak control effort required for the cases with parameter variations is somewhat larger, as expected. It can be reduced by tuning down the control parameters, however, at the expense of roll reduction performance.

Fig. 4 Roll dynamics without parameter variations under dif- ferent controllers: nonlinear feedback (solid), linear feedback (dash—dot), and open loop (dashed)

Fig. 1 The proposed ER damper: (a) Schematic configuration; (b) photograph Here E is the electric field. The a and @ are intrinsic values of the ER fluid to be experimentally determined. In this study, for the ER fluid, the commercial one (Rheobay, TP Al 3565) is used and its yield stress at room temperature is reported by 591E!“? Pa [7]. Here the unit of E is kV/mm. It is herein noted that the Bingham model of the ER fluid [8] is adopted for the derivation of Eqs. (1) and (3).

Fig. 2 Damping force versus piston velocity at various electric fields: (a) measured; (b) simulated

Fig. 4 Bond graph model of the ER suspension system Fig. 3 Mechanical model of the full-car ER suspension system

replacing the sign function by the saturation function with appro- priated boundary layer thickness ¢; for (i= 1,2,3,4). The control- ler u given by Eq. (13) is designed in an active actuating manner. However, the ER damper is a semiactive actuator. Thus, the con- trol input force should be applied to the ER suspension system according to the following actuating condition. In this study, HILS is undertaken to evaluate the performance of the full-car ER suspension system. The configuration of the proposed HILS is schematically shown in Fig. 5. The HILS is divided into three parts: interface, hardware, and software. The interface part is composed of computer (IBM, Pentium-133 MHz) on which DSP board (Mtt, TMS320C31) is mounted. To drive the DSP board in real time, Matlab® program with Simulink® tool- box is used and its looping time is chosen by 0.9 ms. The hard- ware part is composed of the ER damper, high voltage amplifie (Trek 10/10A), and hydraulic damper tester (hydraulic power unit). The software part consists of the theoretical model for the full-car ER suspension system and control algorithm programmed n the computer. The photograph of each component for the HILS is shown in Fig. 6. As a first step, the computer simulation for the full-car ER suspension system is performed with initial value. This computer simulation is incorporated with the hydraulic damper tester which applies the displacement (suspension travel) to the ER damper according to the demand signal obtained from the computer. It is also connected to the high voltage amplifie which applies control electric field determined from control algo- emi

Table 1 Parameters of the full-car ER suspension system rithm to the ER damper. Then, the damping force of the ER damper is measured from the hydraulic damper tester and the measured damping force is fed back into the computer simulation. In short, the computer simultaneously runs both the hydraulic damper tester and high voltage amplifier during simulation loop, and the computer simulation is performed based on the measured data. In this study, the ER damper at rear right side is chosen for the HILS by considering the capacity of the hydraulic damper tester. (=0.856 m/s). The second type of road excitation, normally used to evaluate the frequency response, is a stationary random process [14] with zero mean described by

Fig. 7 Trajectory tracking performance of the hydraulic damper tester: (a) bump and (b) random excitation Here w,=27V/D, z,(=0.035 m) is the half of the bump height, D(=0.8 m) is the width of the bump, D,,,(=2.4 m) is the wheel base which is defined as the distance between the front wheel and the rear one, and V is the vehicle velocity. In the bump excitation, the vehicle travels the bump with constant velocity of 3.08 km/h

Fig. 8 Controlled responses of the ER suspension system for bump excitation: (a) without and (b) with parameter perturba- tions We first investigate the tracking performance of the hydraulic damper tester during simulation loop prior to implementing the HILS. Figure 7 presents suspension travel tracking performance

Fig. 9 Controlled frequency responses of the ER suspension system for random excitation: (a) without and (b) with param- eter perturbations

Let us for simplicity define the time dependent gas velocity of the inlet port, v, (x,,)=v,.. The inlet port gas velocity, v,, can be approximated (Broome [16], Moraal et al. [13]) by the solution to a second order forced differential equation with known initial values:

where, r.=(2/(y+1))”~! is the critical pressure ratio. Mh SSeS ee ee ae ts Fore aR: TSE Fs arses a I SCN eR a ae Sea ae The above expression is simplified by taking into account that for camless unthrottled operation the intake manifold pressure is almost equal to the ambient pressure py= 10° Pa. In addition, we consider 7; =T,=T=316°K, R=287 J/kgk, y=1.4 and that typical operation results in limited backflow through the orifice (p,/p2=1.1). The flow equation

The approximation simplifies significantly the computation an results in fast execution of the simulations. Thus, the flow throug! the throttle body, mg, and the intake valve, m,, are written as where, A,=A,,k (m7) is the effective flow area function for the valve normalized at sonic conditions and constant temperature, and d is defined:

Fig. 1 Camless valvetrain controller (‘inner loop”) and cam- less engine controller (‘“‘outer loop’’)

Fig. 3 Intake valve profiles for different speeds Fig. 2 Electro-hydraulic camless (dashed), simplified camless (solid), and conventional (dotted) intake valve profiles 3 Performance Variables

Fig. 4 Engine breathing characteristics of a 4-cylinder engine at an engine speed of 1500 rpm

Fig. 6 Model (dashed line) versus experimental data (solid line) at 1500 rpm The model described in this paper was evaluated using experi- mental data for a 4-cylinder engine (Moraal et al. [13]) during wide open throttle operation. The valves of the experimental en- gine are cam-driven, therefore, the intake valve profile is the con- ventional sinusoidal profile with the following specifications: (i)

Fig. 5 Cylinder P-V diagrams that demonstrate the pumping loss in a throttled and a camless unthrottled engine (lift profiles are also shown)

Fig. 8 Model (dashed line) versus experimental data (solid line) at 4000 rpm In Figs. 6, 7, and 8 we compare the model with the actual engine data for wide open throttle (f in Eq. (15) is fixed to 90 deg) and engine speeds of 1500, 3000, and 4000 rpm, respec- tively. The manifold and cylinder pressures are plotted for the model (dashed curve) and the experimental engine (solid curve). The model Helmholtz resonator parameters were determined as w,=27*176, €,=0.005 for an engine speed of 1500 rpm, and w,=27*190, €,=0.15 for engine speed of 3000 and 4000 rpm. These values are similar to the ones identified in (Moraal et al. [17]).

Fig. 7 Model (dashed line) versus experimental data (solid line) at 3000 rpm

Note that the model captures most of the dynamic effects. Moreover, the model predicts rather accurately the integral quan- tities, namely, the cylinder air charge and pumping losses, which are overall measures of the engine performance over the induction stroke. The comparison of predicted cylinder air charge (77) and pumping loss (£) versus the actual data cylinder air charge (m) and the pumping loss (£) follows: The maximum error in air charge prediction can be calculated to be less that 10% and the maximum error in pumping losses is 11.3%. Note also, that the model captures the 10 deg overlap between the intake and exhaust valves.

The input-output behavior is derived for fixed intake valve tim- ing, u,, and parametrically varying intake valve maximum lift u, and duration uw, in the crankangle-domain model. We fix u, to zero degrees ATDC, because neither advancing nor retarding u, affects the performance variables in a favorable way (Ashhab et al. [11]). This result is in agreement with the work of Miller et al. [20] where uw, was set to zero at all loads and medium engine speed. Note, that by setting u,=0 deg ATDC, valve closing is approximately equal to the valve duration. The input-output mean-value model can now be described by a static nonlinearity that is identified by fitting data attained by simulation, and a de- lay:

Fig. 11 Investigation of the effects of the higher order inertial and acoustic dynamics to the mean-value model at an engine speed of 1500 rpm Figure 12 shows the discrepancy between the two models for 6000 rpm. At high engine speed, Helmoltz resonator dynamics increase the cylinder air charge during late intake valve closing (ram effect as in Broome [11]). This is achieved with no adverse effect in specific pumping losses. For early valve closing, how- ever, Helmholtz resonator dynamics reduce the cylinder air charge and increase the pumping losses. To explain this behavior we plot part of the P-V diagram for three different intake valve closing timings in Fig. 13. Specifically, Fig. 13 shows the P-V diagrams for early (IVC=80 deg), at the event (IVC=180 deg), and late

(IVC=230 deg) closing timings with and without Helmholtz ef- fects. It is, thus, evident that higher order dynamics in the intake runners introduce considerable uncertainty in the mean-value model representation and have to be considered in the controller development.

Fig. 12 Investigation of the effects of the higher order inertial and acoustic dynamics to the mean-value model at an engine speed of 6000 rpm

Fig. 1 Block diagram of the adaptive feedforward control scheme The feedforward controller is based on the mean-value cylinder breathing characteristics developed in Part I. It calculates the con-

Fig. 2 Cylinder air charge requirements and optimality con- straints for operation at 1500 rpm

We start by summarizing the state-space representation of the multicylinder breathing dynamics. Let the state vector

Fig. 5 Mass air flow and air mass fractions of three-cylinder overlap

Fig. 7 Nominal «< as a function of v, for a four-cylinder en- gine at an engine speed of 6000 rpm Fig. 6 Nominal ,{ as a function of vy for a four-cylinder en- gine at an engine speed of 1500 rmp

where we exploit the fact that m4;43=m4;-1, Ma a3= pai: and pete L— Mii for all k during quasi-static engine condi- tions. Note that Y; is the measurement vector and M; is the cyl- inder air charge vector. To ensure that the periodic conditions are satisfied we perform the estimation only during constant pedal position which does not require changes in valve inputs. Note also that the cylinder air charge vector lags the measurement vector by (r—1) events, where r is the maximum number of cylinders drawing air from the intake manifold (r=2 in four-cylinder case). Estimation of the cylinder air charge vector depends on the values of the entries of matrix Aj, i.e., the values of the mass fractions (u as Values for the yh are calculated eg the nominal simula- tion model. Figures 6 and 7 show nominal jit data as functions of vg and v, at engine speeds of 1500 and 6000 rpm, respectively. For the estimator implementation we ignore the small dependence of uk on vU;. ture of the breathing process allows us to lift the linear time- variant (LTV) system to a 4-input 4-output linear time invariant (LTI) system:

As in the four-cylinder engine case, the observability of the individual cylinder air charge depends on the determinant of the nXn matrix in (22), and consequently, on the values of jk . The following lift and duration values result in conditions for an un- observable system: 4.2 General Case. Individual cylinder air charge estimation for arbitrary number of cylinders, n, and arbitrary (but reasonable r<n/2) number of cylinder overlap, r, is given by

Fig. 8 Flow chart of the closed loop controller algorithm for a four-cylinder engine with two-cylinder overlap (r=2) The closed-loop algorithm for the four-cylinder engine is sum- marized. Consider j=0,1,2..., then where €,= €° if /<0. The conventional engine measurements (iit y and p,,) are used to compute the quantities y,, y., y3, and y4 based on Eq. (12). For simplicity, we use here the equivalent compact state space representation (Eq. (13)): In the next iteration, however, only €, and &, are going to be used. Specifically, €) is not used because cylinder 4 has received already its inputs, w,, and, moreover, y4 has been sampled. In addition, €3 is not updated but kept constant, €;= €_,, in order to satisfy the periodic conditions for having 713=717, and we = ia in the estimation algorithm. The feedforward controller parameters for the next iteration can, thus, be expressed as 5.1 Four-Cylinder Engine Case. With no loss of general- ity, the algorithm is initialized with the input uj=(u jy lt dy) =C(m**, €y) to the fourth cylinder during event 0 (initial event), where &)= €° is calculated based on the nominal engine data (Eq. (6). The first iteration of the controller algorithm (shown in Fig. 8) starts with the feedforward controller calculating the control sig- nals u;=(u,;,Uq.) during events i= 1,2,3,4, that is,

In principle, we can now calculate the controller parameters (Eq. (10): quasy-static periodic conditions in cylinder 4 for which the equali- ties #jp=m,, and pe ey are satisfied, and we can obtain

Fig. 9 Flow chart of the closed-loop controller algorithm for an n-cylinder engine with three-cylinder overlap (r=3)

cycle minus (r—1) fundamental events each iteration of the con- troller algorithm. Thus, two consecutive windows of measure- ments have (r—1) fundamental event overlap.

The simulations are chosen so that we test individually all of the components of the closed loop algorithm, specifically, the adaptive feedforward controller (Section 3) with the cylinder air charge estimation (Section 4). We gradually increase the require- ments introducing model error and cylinder-to-cylinder overlap. Cycle-to-cycle tracking and cylinder-to-cylinder balancing is tested under different scenarios of engine and actuator conditions.

Fig. 10 Nominal engine tracking Fig. 11 Air charge regulation during engine speed changes The control signals (uz and u,) are shown in subplots 2 and 3. The control signals for the four cylinders are identical because the cylinders are balanced. The step changes in the demanded cylin- der air charge forces the feedforward controller to switch between the branches of the ‘‘L’’ shape (see Section 3). In the first engine cycle (t=[0,0.08]), the value of m**’ is larger than the critical air charge (m,,) and thus the feedforward controller selects u, =7 mm (subplot 2) and computes the corresponding u,= 134 deg (subplot 3) that satisfies m“°’. The small difference between the desired and individual cylinder air charge is due to errors in the curve fitting used in the feedforward controller. The closed-loop controller algorithm balances the four cylinders within three en- gine cycles. The next value of m** is less than m,,. Thus, the feedforward controller selects u,z=80 deg (subplot 3) and com- putes u;= 1.36 mm (subplot 2) that satisfies m“*’. As the desired

Fig. 12 Cylinder balancing in the cylinder overlap case

Table 1 Notation in the Moore-Greitzer model

The details are omitted here because of the journal page constraints; see Wang et al. [11] for an outline of the main steps. After extensive calculations, employing MATHEMATICA to solve integrals in closed form,* the compressor model becomes

Fig. 3 Equilibria of the eMG3 model with «=0

Fig. 7 Bifurcation diagrams for the open-loop system with e =0 and £=0.71. The throttle opening y is the bifurcation parameter.

Routine calculations show that in case of e-MG3 one obtains the expressions in Table 2. For a_ general ‘compressor characteristic (not necessarily e-MG3), the critical slopes would be defined as

Fig. 8 Bifurcation diagrams for the open-loop system with e =0.9 and 6=0.71

Fig. 9 Transient responses for throttle opening y=1.15, slightly below the value for the stall inception point. A low value of 6 (@=0.71) results in rotating stall, while a high value of 6 (B=1.6) results in surge.

Let us consider an equilibrium on the stall characteristic. The equilibrium is determined by the value of I’, which is given as The controller development in this section is independent of the form of compressor characteristic. We only require that V (1) =0 and V&(1)<0, ie., that V(®) has a maximum at ®= 1.

Table 3. The stability conditions for the system (24)—(26)

Fig.10 The function \V,(W). Variations in I result in no more than one stall and one axisymmetric equilibrium. A qualitative sketch of My(‘V) is given in Fig. 10. This function is chosen as monotonically nonincreasing. In presence of varia- tions in I, MVy(‘V) will ensure that there is never more than one stall equilibrium and one axisymmetric equilibrium. If My(V) were linear in V, this would allow a creation of a second axisym- metric equilibrium.

For a cubic characteristic, these conditions have already been re- vealed by Krener [14]. Since these conditions involve not only the shape parameters S$, and S, of the compressor characteristic but also the value of pressure at the peak, 2+W eo, we are of the opinion that controller (23), in which I’ is the bifurcation param- eter, is preferable to the controller (98). where yo is a set point/disturbance parameter and the gains kp, kp, ky, and kq are required to satisfy For right-skew compressors, for which S,;>0 and S,<0, the conditions (99)—(101), in particular, allow controllers of the fol- lowing types:

where &(R)=B?/7VVp,(R). The characteristic polynomial of the Jacobian (117) is Sufficient conditions for stability of the system near the stall in- ception point are readily deduced from (118) and are given in Table 4. A center manifold argument similar to that in Section 6.3 proves stability at the bifurcation point.

8.2 Actuator/Sensor Trade-Off. The time constant 7 ap- pears only in the third condition in Table 4. When cg=0, this condition becomes We point out that in the case of a linear implementation (98) the gains kr, kg», kw, and k@ are required to satisfy condition: (99) and (101), whereas the condition (100) is modified as

Fig. 13 Bifurcation diagrams for the case e=0.9, B=0.71, with the (W,R)-controller. The control gains are cp=3 and cy=1.

Fig. 14 Bifurcation diagrams for the case «e=0.9, G=0.71, with the nonlinear (W,@)-controller. The controls are gain dy=2 and nonlinear function Vy(W).

Fig. 15 Bifurcation diagrams for the case e=0.9, B=1.42, 7=0.44, with the full-state controller. The control gains are Crp=43, Cy=17, and Cy=22.

Table 5 h,(0) for nonlinear implementation Thus, high gain on ® makes a right-skew compressor unstable. Similar but more complicated argument applies to h4(Ro).

Table 7 h,(0) for linear implementation Table 8 h; for linear implementation A discussion similar to that in Section 10.1 shows that the increase of kq results in instability for the right-skew case. Let us now set kp=ky=0 and concentrate on the UTRC controller. Since h,(0) and Ry are independent of kg, for stability we re- quire that h5(0)>0. It is obvious that, for large kg, the sign of h3(0) is determined by the sign of S,S)=—4¥U(1). Thus, the

Since Ro is independent of dg, the sign of h,(0)+h;(0)Ro for large dp is determined by the sign of S,;S,= —4 c(1). This means that, for large do , - For small Ro>O (near the stall inception point), a necessary condition for stability is h,(0)+h;(0)Rj>O. Let us concentrate on The quantities 4,(0) and 4;(0) are given in Tables 5 and 6, re- spectively.

Fig. 16 Closed-loop bifurcation diagrams for e=0.6 with UTRC controller. The control gains are kKp=3 and kg =1,10,100.

The pressure Py can be obtained by substituting the solution (36) into the Eq. (1b), which gives Similarly, the approximate solution of problem (18) is the same as the approximate solution of problem (10), that is, the transfer functions are given by Eqs. (27)—(28). Finally, we consider the approximate solution of problem (19), where flow rate is given at beginning of the line and pressure is given at end of the line. The flow rate solution is given by

Table 1 Some window functions If the trigonometric basis (32) is used with the function

Fig. 1 Imaginary part of the propagation operator ['2(s =ifwT) and its approximations, when friction coefficient « =1/10

The Q model is used when transmission lines are connected together or to fluid volumes (e.g., hydraulic cylinder) by compo- nents of negligible volume such as directional control valves (Fig. 2). Such connectors are modeled as orifices as in Piché & Ellman [18] and Ellman and Piché [19]. The Q-model transfer functions are A 4-mode SIMULINK realization of the Q model is shown in Fig. 3.” Implementation in other simulation environments should be straight forward. Because the models are decoupled, the model has an obvious parallellization. The approximate transfer function (59) are rational polynomials thus inverse Laplace transformation can be calculated giving constant coefficient ordinary differential equations. For ODE-based simulators the following state space realization can be used (with n modes):

Fig. 3 SIMULINK realization of the Q-model with 4 modes 2SIMULINK realization of Q, P, and PQ models of transmission lines are avail- able at ftp.cc.tut.fi/pub/math/piche/fluidpower

Fig.2 A pipe system with Q-pipe models. Pipes are connected by the orifice that simply calculates flow rate through it when pressure drop is known.

Table 2. The physical properties of pipe line system

The PQ model is used when the connections at the ends of the transmission line are of different types. One end can be connected to an orifice and the other directly to a volume. The PO-model transfer functions are The P model is used when transmission lines are connected together to fluid volumes by components of negligible resistance (Fig. 4). Such connectors are modeled as hydraulic volumes (Piché and Ellman [7]). The P-model transfer functions are

Fig. 4 A pipe system with P-pipe models. Pipes are connected to the volume that is modeled by an integrator. More than two pipes may be connected to a single volume.

Pressure in time domain at intermediate points along the transmis- sion line can be calculated using the Ritz approximation (23)

When modeling a network containing several transmission line elements, one should select the number of modes such that the ratio n/L is about the same for each element. This ratio reflects the bandwidth over which the model is supposed to be accurate. Thus, longer transmission lines require more modes to accurately model the same frequency range as shorter lines with fewer modes. Assigning too many modes to short transmission lines is inefficient, since the numerical ODE integrator will need to use small time steps to follow the high frequency modes. This numeri- cal cost would be especially evident with simulators that use BDF methods (such as Gear’s method) as the integrator, because of their inefficiency with oscillatory problems.

Fig. 6 Single pipe example, comparing dissipative (2D vis- cous) and linear friction (1D viscous) models of pressure re- sponse at beginning of line

Fig.5 Single pipe example, pressure response at beginning of line

Fig. 7 Water hammer example pressure response at valve for first cycle of water hammer in high-viscosity oil Table 4 The physical properties of water hammer system in Sl units (Holmboe and Rouleau [20])

Substituting the Ritz approximation into Eq. (70) and simplifying yields the ordinary differential equation system of second order Equations (66a) and (66b) do not include frequency dependent friction effect, since the equations are one-dimensional. Equations (71)—(73) can be modified to include frequency dependent friction effects (2D friction), if we assume that for small flow rate and pressure perturbations the turbulence flow behave the same way as laminar flow. This assumption is related to turbulent mean flow condition in Trikha, [12]. Comparing the Q model (59) and Eq. (72) and (73) we can modify the matrices as follows (2D friction):

Fig. 8 Nonlinear (turbulent) flow example pressure response at beginning of line

Fig. 1 Diagram of a pneumatic vibration exciter It was experimentally found that in real flow conditions present in a working generator temperature changes can be omitted and Pressure p,, controlled by valves f; and f,,, mainly depends on the quantity of air mass in chamber V,. Changes in the chamber volume caused by piston vibrations have little effect on pressure P1, in the open chamber and air stiffness in such a chamber is, inter alia, dependent on frequency w and is smaller than in the closed chamber [1]. The component of pressure function pj, re- sulting from varying volume of chamber V, is a feedback signal in the generator. When the feedback is compensated, changes in the chamber volume do not cause the air pressure in this chamber to change, thus the air temperature remains unchanged, too. This leads to a conclusion that the air temperature 7, in chamber V, depends mainly (or exclusively after feedback compensation) on air flow processes. _ «2 a el oe In ee

Fig. 2 Physical model of a pneumatic vibration exciter

Assuming the initial conditions g,(0)=0, q2(0)=0, and q3(0) =0, the Eqs. (24) can be solved numerically by one of the The system of Eq. (23) can be written in the state coordinates assuming the notations x=q,, X=q2o, X=q2o, p=q;3; hence

into a variable component chamber /0 and constant component chamber //. There is a cylinder linear in the variable component chamber. The liner consists of four ring segments /2, 13, 14, and 15. Between them are rectangular slots: an inlet 76 and outlet 17. These slots have rectangular cross-sections. Their longer sides with length / are parallel to the cylinder axis. The cylinder liner design allows each segment to turn independently. This, in turn, enables separate width adjustment of both slots b, and b, and the phase angle \, between them (Fig. 5). In Fig. 4 and 5 only one inlet and one outlet slots are shown. In fact, there may be more such slots as described in the earlier work; and they are arranged symmetrically on the cylinder liner wall. If this is the case, the opening functions f; and f,,,; are formed from a number of smaller components. This facilitates achieving the desired character of air flows and enhances the effectiveness of generator operation. The rotor /8, placed in the cylinder liner, is used for harmonic control of the opening function. Its head surfaces have sinusoidal contours with the amplitude A, and nonfractional and identical number of periods on each side. The phase shift angle of the sinusoidal con- tours is zero. The rotor is turned, by a device not shown in the drawing, at a variable (rotating) speed n dependent on the required exciting frequency w. The cylinder liner design with sinusoidal contours of the rotor is shown in Fig. 5. The values of the variable components of the opening functions depend on the angular posi- Fig. 4 Pneumatic vibration generator design for experimental investigation. 7, cylinder; 2,3, covers; 4, two-sided piston; 6,8, inlet holes; 7,9, outlet holes; 10,17, chambers; 12,13, 14,15, ring components of cylinder liner; 76, inlet slot; 77, outlet slot; 78, rotor with sinusoidal contours. Fig.5 A fragment of the extended cylinder liner with inlet and outlet slots, cylindrical rotor surface with sinusoidal contours marked with a broken line; 1,2,3,4, ring components of cylinder liner; 5, inlet slot; 6, outlet slot; 7, rotor with sinusoidal contour

Fig. 3 Variable of a pneumatic exciter state

Fig. 4 Pressure signal versus time for G=0, V,g=30 mm’, F=5Hz; (a) L,=150 wm; V=105 mm’; (b) L,=150 wm V=350 mm; (c) L,=150 wm; V=665 mm; (d) L,=45 wm; V=675 mm*; — experimental data; - - theoretical data

Fig. 2 Geometrical configuration of an OPD with a grid system

Fig. 5 Pressure signal versus time for G=0.5*10-° m%/s Pa, Vag=70mm*: (a) L,=150 um; V=105mm'; F=5 Hz; (b) L,=150 wm; V=200 mm’; F=5 Hz; (c) L,=150 wm; V=350 mm’; F=10 Hz; (d) L,=150 wm; V=665 mm; F=10Hz; — experi- mental data; - - theoretical data

where Ax,At are step sizes in the space and time domain, respec- tively, superscript ‘‘’’ stands for the time step. From Eqs. (4), (7), and (8) we obtain the equation for defining the temperature at the exposed absorbent surface at each time (t+At) step

Fig. 7 Pressure signal versus time for G=0.5*10-° m/s Pa, Vag= 70 mm, L,=150 wm; V=105 mm®: (a) F=10 Hz; (b) F=50 Hz; (c) F=100 Hz; (d) F=200 Hz; — experimental data; - - theoretical data

Fig. 1 Experimental apparatus for the temperature measure- ment (during charging)

Fig. 2 Experimental apparatus for the temperature measure- ment (during discharging)

Fig. 3 Temperature responses during discharge

Fig. 4 Temperature responses during discharge (a) tempera- ture response assuming that the whole heat generated by the expansion of air use transferred to the material

Fig. 6 Experimental results of the steady flow rate measure- ment 5.2. Results and Discussion. Figure 6 shows the pressure curves and the flow rates of the experimental results. The experi- ment was done for four flow rates. Cases 1 to 4 correspond to supply pressures of 542, 444, 346, and 248 kPa, respectively. The

Fig. 7 Experimental apparatus for the unsteady flow rate mea- surement

Fig. 8 Unsteady flow rate measurement with the isothermal chamber (f=5 [Hz])

Fig.9 Unsteady flow rate measurement with isothermal cham- ber (f=40 [Hz])

Fig. 10 Unsteady flow rate measurement with the normal chamber (f= 40 [Hz])

The first term in Eq. (2) is the driving moment from the torque- motor. The last five terms are restoring moments on the armature from the net rotational stiffness, the cantilever feedback spring stiffness (if this spring exists), the damping in the torque-motor, the pressure difference across the nozzle, and the dynamic flow force on the flapper, respectively.

The first two terms in Eq. (6) are the forces from the spool pres- sures and the spool damping, respectively. The next two terms in Eq. (6) are, respectively, the flow reaction force, which depends on the actuator pressures, and the transient flow force, which de- pends on the actuator flow rates. These two terms show the de- pendence of the servovalve on the actuator dynamics given below. The last term in Eq. (6) is the restoring force from the feedback spring, where

Fig. 3 Frequency responses of the servovalve The experimental system consists of a Moog servovalve utiliz- ing a cantilever feedback spring and a 0.9 kg double-ended actua- tor with an effective area of 3.613 e * m? (0.56 in’). The supply pressure is 18.6 MPa (2700 psi). Mobil DTE Light oil is used at an average temperature of 90°F. The actuator and spool are con- nected to linear variable differential transducers (LVDT) with +3 mV noise level. The range of the actuator LVDT is +2.5 V for the total stroke of 0.0254 m (+1 in). While the servo valve physical parameters are not available for model verification, the experi- mental frequency responses for the servovalve and the actuator are measured with a signal analyzer, using a ‘“‘swept sine’’ method that generates fixed-amplitude sine waves of varying fre- quencies. This is then used to validate the model order and struc- ture. It is desired to measure this data near the equilibrium actua- tor position x,,=0. Since the actuator is of type 1 as Eq. (22) indicates, it is difficult to keep the actuator near the equilibrium position in open loop. Therefore, a proportional control loop (gain= 1) is used to create a stabilized plant, as shown in Fig. 2. With different input amplitude levels (15 mV, 30 mV, and 50 mV), the frequency responses of the servovalve and actuator are measured and shown in Figs. 3 and 4, respectively. Using the frequency responses from the three input amplitudes, an averaged frequency response is computed for the servovalve and the actua- tor, and nominal models are fitted through the use of an equation- error method [17].

Fig. 4 Frequency responses of the actuator (dot line: 15 m\ input; dash line: 50 mV input; solid line: model) The nominal model for the averaged actuator frequency re- sponse is found to be

Fig. 5 Comparison of servovalve models (solid line: average of experimental data; dot line: second order model; dash line: third order model; dot-dash line: derived model) not vary much with respect to input amplitude. The averaged ser- vovalve frequency response is shown in Fig. 5 with the best curve fits for the models of Eqs. (24)—(26). The derived servovalve model of Eq. (24) is

Fig. 6 Mixed sensitivity problem as a standard problem The solution of the mixed sensitivity problem provides an ap- proximate solution of the robust performance problem because for any complex numbers ‘‘a’’ and “‘b,”’

Fig. 7 Multiplicative uncertainties of the plant models (dash- dot line: W,;; dot line: W,¢; dash line: W,,)

Fig. 9 Experimental step responses (dash-dot line: using Ks; dot line: using K,; dash line: using K,; solid line: plant) Figure 9 shows the experimental step response using each ro- bust performance controller for a reference amplitude of 30 mV. Better system response is achieved with the more accurate eighth-

Fig.8 Performance weights and experimental sensitivites (dot line: using K, from the second order model; dash line: using K, from the third order model; solid line: using Kg from the derived fifth order model)

Recall that K,=0 for a servovalve with a cantilever feedback spring and K=O for a servovalve with a direct spool feedback spring. ~ The three equations ((Al), (A2), and (A3)) contain three un- knowns: p,)(5)—Ps2(s), O(s), and x,(s). These three unknowns are solved as a function of the input current.

Fig. 1 Adaptive control system block diagram Then the plant equation in (6) becomes

Adaptive Repetitive Control. Since both the reference tra- ectory and disturbance are periodic, S(q~')=1—q~™ may be Note that some coefficients are comparably smaller than others. It is possible to reduce the nominal model to a lower order one. However, to achieve closed-loop stability and better tracking per- formance, we use the 8th-order nominal model so that the size of the unmodeled dynamics « is sufficiently small as suggested by Theorem 1. Indeed, the adaptive control experiment was unstable when using lower order model for the same tracking reference profile.

A cam profile shown in Fig. 3 is used as the tracking reference signal throughout the experiments. To achieve high accuracy, the cam profile consists of 256 points over one spindle rotation, ren- dering 2560 Hz sampling frequency at 600 rpm spindle speed. Throughout the experiment the fluid supply pressure is 20684.4 Kpa (3000 psi) unless otherwise mentioned. Off-line System Identification. Five sets of frequency re- sponse were obtained for the actuator with different actuation am- plitudes by the swept sine method (Fig. 4). The curve fit model shown as solid line in Fig. 4 shows unmodeled dynamics from the real system. The curve fit discrete model considered as the tuned nominal model in Eq. (2) is

Fig. 4 Swept-sine frequency responses at different input magnitudes

Fig. 5 Tracking error and control signal for nonadaptive re- petitive control with perturbed nominal model

Fig. 7 Steady-state tracking error and control signal for adap- tive repetitive control with perturbed nominal model Fig.6 Transient tracking error and control signal for adaptive repetitive control with perturbed nominal model

Fig. 8 Tracking errors of nonadaptive and adaptive repetitive control at 13789.6 Kpa (2000 psi)

Fig. 9 Tracking errors of nonadaptive and adaptive repetitive control at 10342.2 Kpa (1500 psi) The order of the closed-loop characteristic polynomial for the repetitive control system is usually high and has poles close to the unit circle. Thus the closed-loop system becomes unstable with the slight perturbation. All industrial actuators have manufactur- ing variations. If a fixed repetitive controller is applied to several actuators, it may not be stable due to parameter mismatch. An adaptive repetitive controller can handle the parameter mismatch and adapt itself to new controller gains.

Since the right-hand side of the above inequality is nondecreas- ing, we have

(II): r(k)=0 since S(q~')y,e(k)=0. Also since ||¢(k)|| is bounded from (I), we have D(k) €1,5D €1,>D(k)—0. With X=max(A.Ay), (A2) and (A3) imply that The RHS of the above inequality is nondecreasing. Provided that e<X(1—X)/(k,||v|]) we have Further, § has no common factor with B, by assumption, there exist polynomials 7 and A, such that yS+AB,=1. Therefore, xSu+ABu=u and xi+A(A,y—C,d)=u. Therefore u(k) is bounded since the left hand side of the above is bounded. Since || @(n)||Skgl|x(n)|| and D(n) €1,, by the discrete form of the Bellman-Gronwall’s lemma (see Lemma 3), ||x(k)|| is bounded. Therefore # and y are bounded.

Under ideal conditions, where the load torque is zero and sys- tem parameters such as inertia load and supply pressure are con- stant, the required response may be achieved by the proper design of a controller using the pole-placement technique. Hence, con- ventional controllers are designed to such nominal conditions. However, in real systems, load torque disturbances and variations in system parameters, such as inertia load and supply pressure, The transfer functions of the speed response for input command and external load torque are obtained as follows from Fig. 3: Figure 3 shows a block diagram of the system with a conven- tional cascade speed controller. The basic control structure is a double loop, which consists of an inner proportional position- control loop and an outer proportional and integral speed-control loop.

Fig. 3 Block diagram of the conventional control system

Fig. 1 Schematic diagram of VDHM control system

Fig. 4 Block diagram of the proposed control system

Table 1 System parameters Figure 7 is a schematic diagram illustrating the experimental setup. It consists of three parts: a motor system, a load generating system, and a data acquisition and control system. Some of the technical specifications of the experimental setup are as follows: the sampling period is 1 ms, the servovalve rated current is 60 mA, and the rated flow is 30 l/min at a pressure drop of 70 bar. The max. displacement volume of the VDHM is 40 cc/rev, the encoder resolution and the A/D, D/A converters are 1000 pulse/ rev and 12 bits, respectively. Bandwidth of the amplifier for the LVDT is 500 Hz and that of the signal conditioner for the torque sensor is 15 kHz. The signal of the motor speed from the encoder is so noisy that a low-pass filter with a bandwidth of 50 Hz is used. The nominal operating conditions of input speed command, inertia load, and supply pressure are 1000 rpm, 0.04 kgm? (8 times of motor shaft inertia J,), and 150 bar, respectively. Some of the experimental conditions are given as follows:

The transfer function G(s) can be obtained from Fig. 4 as fol lows:

Fig. 8 Comparison of measured and observed torque

Fig. 10 Response for the sinusoidal load torque Fig. 9 Response for the step load torque

Table 2 Overshoot for various 7, and 7, (J=2J,).

Fig. 11 Step response for the input change

Fig. 12 Step response for the change of inertia (dotted line: conventional PI controller, real line: proposed controller)

It is herein noted that the yield stress of the ER fluid employed in this work gradually increases as the temperature increases to about 90°C, and sharply decreases above 90°C, as shown in Fig. 1(b). This is due to that as the temperature increases, the polar- ization effect of the conductive particles becomes more active. This microscopic phenomenon results in the increment of the yield stress. However, it is observed that the field-dependent ER effect significantly degrades above 90°C, since the particles of the employed ER fluid start melting around 90°C. The information on the variation of the ER effect with respect to the temperature is

Fig. 1 Bingham property of the employed ER fluid. (a) Shear stress, (b) yield stress.

Fig. 3 The manufactured cylindrical ER valve. (a) Assembly drawing, (b) photograph.

Fig. 4 Field-dependent step responses of the ER valve

Now, from the above equations, a state-space model for the governing equation of motion of the proposed ER valve bridge- cylinder system can be written as follows:

Fig. 7 A neural network controller incorporating the PID con- troller. (a) Block-diagram, (b) network architecture.

Fig. 9 Set position control responses without cart. (a) Simu- lated, (b) measured. component of the cart is numerically obtained by differentiating the displacement signal. The pump makes a flow of the ER fluid, and thus the piston rod of the cylinder is moved. The temperature of the operating ER fluid is maintained at 55°C during control action by operating the air cooling system. The flow rate of the ER fluid from the pump is chosen by 8 I/min. The displacement of the cart is measured and compared with the desired displacement. In order to minimize the position error, the control electric field

Fig. 8 Experimental configuration for the position control. (a) Schematic diagram, (b) photograph.

Fig. 10 Set position control responses with cart. (a) Simu- lated, (b) measured.

Fig. 11 Square trajectory tracking control responses with cart. (a) Simulated, (b) measured.

Fig. 12 Sinusoidal trajectory tracking control responses with cart. (a) Simulated, (b) measured.

Fig. 14 Measured tracking durability of the control system Fig. 13 Comparison of measured tracking responses with cart. (a) Neural network, (b) PID.

Fig. 1 Cylinder head arrangement of a digital-displacement pump-motor Figure 3 shows a schematic representation of a motoring cycle [4]. To enable a cylinder for motoring the controller closes the low-pressure valve shortly before the piston reaches the top-dead- center (TDC). Once the valve is closed, the cylinder pressure rises to equal to that of the high-pressure manifold by the time it reaches TDC. The high-pressure valve can then be opened and latched. The piston is then propelled by the fluid pressure toward BDC. In a similar fashion to the TDC valve-sequence, the high- pressure valve is closed prior to BDC such that the residual piston stroke can de-pressurize the cylinder and allow the low-pressure

Fig. 2 Change of the enabling sequences of a 10 cylinder pump and its effect on flow as the demand drops to 3/4th of the capacity

Fig. 4 Pump response at 10 percent to 90 percent step rise of flow demand at a constant speed (pressure-control mode) The displacement out of the accumulator S$, can be calculated from the pressure difference between the current and the previous decisions:

Fig. 5 Pump response to 17 percent speed change, while the flow demand remains constant (pressure-control mode)

Fig.6 Pump-motor following a variable demand curve at 1800 rpm (flow-control) Pave=a running average of pressure of last 5 decisions

The digital-displacement technology has advanced in the last few years both in terms of compu ter simulation and the building of prototypes. The first prototype of a six-cylinder pump was built in 1990. The single cylinder pump-motor prototype (1994) was followed by a six-cylinder axial pump-motor (1996). Figure 13 shows a photograph of the various components used for a six- cylinder radial pump-motor whic writing. Experimental work has simulations, some of which have Figure 14 shows the experimenta motor, running at 1500 rpm und h is under test at the time of shown good agreement with already been published [3,4]. results of a 6-cylinder pump- er pressure-control mode. The pressure is held constant (about 70 bars) while the pump flow follows an increase from 10 to 90 percent of the capacity in about 30 ms. The experimentally measured response time of the digital- Fig. 9 Power requirement of the pump-motor bank

ry [he second bank also consists of a 10 cylinder machine, which only acts as a pump and which is sized to provide the maximum flow demanded by the load. The cylinder enabling sequence is dictated to follow the flow requirement of the load, converting the shaft output to hydraulic power. Figure 11 shows the output flow from the second bank as well as the flow-demand. Figure 12 shows the power output from the system compared to the load- demand. It also shows the much smaller constant power contribu- tion of the prime-mover.

Fig. 13 Components of a 6-cylinder pump-motor

Fig. 11 Demand and flow from second bank (pump)

Fig. 1 General machine configuration While the valve plate is held in a fixed position, the coupled shaft and cylinder block are driven about the x-axis at a angular speed, w, which will also be considered a constant in the follow- ing analysis. During this motion, each piston periodically passes over the discharge and intake ports on the valve plate. Further- more, because the slippers are held against the inclined plane of the swash plate, the pistons undergo an oscillatory displacement in and out of the cylinder block. As the pistons pass over the intake port, the piston withdraws from the cylinder block and fluid is drawn into the piston bore. As the pistons pass over the discharge port, the piston advances into the cylinder block and fluid is

where A, is the effective pressurized area within a single piston bore of the cylinder block and will be discussed later. Similar to Eq. (1), the forces acting on the cylinder block in the y- and z-directions may be summed and set equal to zero as follows:

where E,, and E, locate the centroid of the valve-plate reaction force, F,,. Summing the moments acting on the cylinder block about the origin of the y- and z-axis and setting them equal to zero yields the following results: Piston-Slipper Assembly Free-Body-Diagram. The forces acting on a single piston-slipper assembly are shown in Fig. 3. These forces result from the pressure within the nth piston bore (A,P,,), the equal and opposite reaction of the cylinder block against the nth piston (—F), and —F7, ), and the net reaction between the nth slipper and the swash plate (F sw, )* Summing the forces acting on the nth piston-slipper assembly in the x-direction and setting them equal to the time rate-of-change of linear mo- mentum for the entire assembly in this direction yields the follow- ing result:

Fig. 3 Piston-slipper assembly free-body-diagram

Fig. 4 Schematic of a piston-bore control volume For a constant speed machine dt=d6/w. Using this result with Eqs. (12), (13), and (14), the pressure rise-rate within a single piston chamber may be written as

Equation (15) is a nonlinear, first-order, differential equation that must be solved numerically. Figure 5 shows a typical numeri- cal result of Eq. (15). Note that both pressure, P, and port area, A,, are plotted in this figure. As shown in Fig. 5, the typical numerical solution to Eq. (15) demonstrates rather uninteresting behavior for the pressure, P. As the piston bore passes over either the intake port or the discharge port of the valve plate, the port area, A, , remains at a maximum constant. Within these regions, the pressure within the nth piston-bore also appears to remain fairly constant (i.e, P=P; or Pz). The two ports on the valve plate are bridged by transition regions where A, goes from a maximum value to a minimum value, slowly grows within the transition slot, and then quickly returns to the original maximum value. As the nth piston-bore passes over the transition regions, the pressure changes almost linearly from one port pressure to the other. Fig. 5 A numerical pressure-profile for the nth piston exhibit- ing essentially no overshoot or undershoot in the transition regions

Fig.6 A numerical pressure-profile for the nth piston exhibit- ing overshoot and undershoot in the transition regions Figure 6 shows another result of this study where the pressure drop between ports has been reduced. Figure 6 represents a run using one sixth of the discharge pressure of Fig. 5. From Fig. 6 it can be seen that a lower pressure drop between ports tends to create significant pressure spikes within the transition regions of the valve plate. The reason for this peculiarity is strictly a result of the volumetric compression and expansion within the chamber. In the first case, the chamber volume decreases at a rate faster than the fluid can squeeze out through the port. If the boundary pres- sure is not sufficiently large compared to the starting pressure, volumetric compression of the fluid will cause the pressure within the piston bore to overshoot the approaching boundary condition.