# Nonlinear Mathematical Modeling in Pneumatic Servo Position Applications

Authors: Antonio Carlos Valdiero; Carla Silvane Ritter; Cláudio Fernando Rios; Marat Rafikov
Journal: Mathematical Problems in Engineering
 (2011)

---

## Abstract

This paper addresses a new methodology for servo pneumatic actuators mathematical modeling and selection from the dynamic behavior study in engineering applications. The pneumatic actuator is very common in industrial application because it has the following advantages: its maintenance is easy and simple, with relatively low cost, self-cooling properties, good power density (power/dimension rate), fast acting with high accelerations, and installation flexibility. The proposed fifth-order nonlinear mathematical model represents the main characteristics of this nonlinear dynamic system, as servo valve dead zone, air flow-pressure relationship through valve orifice, air compressibility, and friction effects between contact surfaces in actuator seals. Simulation results show the dynamic performance for different pneumatic cylinders in order to see which features contribute to a better behavior of the system. The knowledge of this behavior allows an appropriate choice of pneumatic actuator, mainly contributing to the success of their precise control in several applications.


---

## Body

## 1. Introduction

This work presents a new methodology to identify the main nonlinear characteristics in pneumatic actuators and its mathematical modeling in engineering applications. The pneumatic actuator is very common in industrial application [1] because it has the following advantages: its maintenance is easy and simple, with relatively low cost, self-cooling properties, good power density (power/dimension rate), fast acting with high accelerations [2] and installation flexibility. Also, compressed air is available in almost all industrial plants [3].However, there are difficulties of control due to various nonlinear characteristics of the system [4, 5]. The nonlinearities present in pneumatic actuators are motivated by its very low stiffness (caused by air compressibility), inherently nonlinear behavior, parameter variations, and low damping of the actuator systems, which make it difficult to achieve precise motion control. The main nonlinearities in pneumatic servo systems are the servo valve dead zone [6], air flow-pressure relationship through valve orifice [1, 7], the air compressibility, and friction effects between contact surfaces in actuator seals [8, 9].Several recent authors present a study on the characteristics of nonlinear pneumatic actuators [1, 5, 7, 9–11]. Valdiero et al. [6] present a mathematical model to dead zone in pneumatic servo valves, followed by the method used for compensation that is made with the addition of an inverse dead zone function in the control system. Rao and Bone [1] present a modeling approach where they use the bipolynomial functions to model the valve flow rates, but a poor classical friction model is used. Perondi [10] developed a nonlinear accurate model of a pneumatic servo drive with friction, where the nonlinear airflow relationship between the pneumatic valve’s driving voltage and the upstream/downstream pressures is proposed. Endler [7] used the methodology of optimal feedback control for nonlinear systems proposed by Rafikov et al. [12] in servo pneumatic system, and simulation results show that a full nonlinear mathematical model is important in pneumatic robot applications.The main paper contribution is to systematize its nonlinear mathematical model with some innovations such as a new equation for valve flow rate and to show how it is important for the success in control applications. The paper is organized as follows. Section2 brings a description of servo pneumatic positioning system with its main components, the used test rig, and a schematic drawing with the nonlinearities present in the actuator. In Section 3, the systematic methodology of the pneumatic actuator nonlinear mathematical modeling is shown. Results are presented in Section 4. Conclusions are outlined in Section 5.


## 2. Pneumatic Servo Position System

The servo pneumatic positioning system considered in this paper is formed by a proportional servo valve (component 4 in Figure1) and a double-action rodless cylinder (component 2 in Figure 1). This actuator permits to position one load in a desired position of the actuator curse or follow a desired trajectory. Figure 1 shows the schematic drawing of used experimental setup with main components for the purpose of investigating the nonlinear mathematic model. The used acquisition and control system is a dSPACE DS 1102 board. It is composed by 4 analog inputs (ADCs) and 4 analog outputs (DACs) as shown in Dspace [13]. Sensors permit to measure air system inlet pressure (1), the actuator position (3), and actuator chamber pressures (pa and pb) (5) and (6).Figure 1
Experimental setup with main components.Figure2 shows the schematic drawing of a servo pneumatic actuator for better understanding of system behavior. During the operation, the control signal u energizes valve’s solenoid, so that a resulting magnetic force is applied in the valve’s spool, producing the spool displacement. The spool displacement opens control orifices, so that one port is connected to the supply’s pressure line, and the other is connected to the atmosphere. Consequently, there is the pressure difference between cylinder chambers, resulting in a force that moves the mass M in a positive or negative displacement y, depending on the control signal input.Figure 2
Schematic drawing of a pneumatic servo system.Figure3 shows block diagram of the main dynamics in the nonlinear mathematical model of the pneumatic actuator. The main nonlinear characteristics of this dynamic system are servo valve dead zone, air flow-pressure relationship through valve orifice, air compressibility, and friction effects between contact surfaces in actuator seals.Figure 3
Block diagram of the main nonlinear dynamics in the mathematical model.Dead zone is common in pneumatic valves because the spool blocks the valve orifices with some overlap, so that for a range of spool positions, there is no air flow [6]. It is located at the dynamic system as a block diagram shown in Figure 3 and is characterized in Section 3.1.The air flow-pressure relationship through valve orifice is a nonlinear function that depends on pressure difference across the valve orifice and valve opening [7]. In this paper, we present a new mass flow rate equation in Section 3.2.The pressures dynamic model is obtained from continuity equation and results in nonlinear first-order differential equation. This dynamic behavior depends on pneumatic cylinder size. Small cylinder bore size produces significant effects such as a faster pressure response [1]. If the bore size is reduced, the ratio of friction force to maximum pneumatic force increases, and the chamber pressures are more sensitive to small variations in the mass flow rate. Therefore, the precise tracking control is more difficult with smaller bore sizes. This detailed nonlinear dynamics is presented in Section 3.3.The nonlinear friction is the most important factor that affects the motion equation. Friction is a nonlinear phenomenon which is difficult to describe analytically [8]. The friction often changes with time and may depend on an unknown way of environmental factors, such as temperature and lubricant condition. Even so, the modeling of their main characteristics is important. In this paper, we consider the actuator friction dynamics described by the LuGre model, proposed in Canudas de Wit et al. [14], and improved by Dupont et al. [15] in order to include stiction effects. This model is presented in Section 3.4.


## 3. Nonlinear Mathematical Modeling

The systematic methodology of the pneumatic actuator nonlinear mathematical modeling is presented from experimental data and recent literature information. The full system constitutes a fifth-order nonlinear dynamic model of the pneumatic positioning system and considers the nonlinearity of the dead zone, the mass flow rate, the pressure dynamics, and the motion equation, that includes the friction dynamics.


### 3.1. Dead Zone Nonlinearity

This section presents the mathematical model for dead zone nonlinearity and its graphical representation. Dead zone is a static input-output relationship which for a range of input values gives no output. Figure4 shows a sectional view sketch of typical spool valve with main mechanical elements.Figure 4
Sectional view sketch of typical spool valve with main mechanical elements of the proportional valve with input dead zone.The mathematical model for dead zone in pneumatic servo valves presented in this section was obtained from Tao and Kokotović [16]. The dead zone analytical expression is given by the equation(3.1)uzm(t)={md(u(t)-zmd)ifu(t)≥zmd,0ifzme<u(t)<zmd,me(u(t)-zme)ifu(t)≤zme,
where u is the input value, uzm is the output value, zmd is the right limit of the dead zone, zme is the left limit of the dead zone, md is the right slope of the output, and me is the left slope of the output.Figure5 shows a graphical representation of the dead zone. In general, neither the break-points (zmd and zme) nor the slopes (md and me) are equal.Figure 5
Graphical representation of the dead zone.In current fluid power literature, dead zone in valves is expressed as a percentual of spool displacement. Valdiero et al. [6] present a new methodology for dead zone nonlinearity identification in proportional directional pneumatic valves. It is based on observing the dynamic behavior of the pressure in the valve gaps. The dead zone nonlinearity is a key factor that limits both static and dynamic performance in feedback control of fluid power systems. The usual method to cancel the harmful effects of dead zone is to add its fixed inverse function into the controller. This inverse is modeled by a set of parameters that need to be identified. The classic dead zone parameter identification uses expensive flow transducers and special test rig, while our proposed methodology needs only the pressure transducers shown in Figure 1. Experimental results are presented in Valdiero et al. [6] and illustrate the efficacy of this methodology that is cheaper and faster.


### 3.2. Mass Flow Rate

According to Rao and Bone [1], the mass flow rate model of the proportional valve is a key part of the system model. In this paper, we use an innovator model to mass flow rate equation qma and qmb developed by Endler [7], given by equations(3.2)qma(u,pa)=g1(pa,sign⁡(u))arctg(2u),qmb(u,pb)=g2(pb,sign⁡(u))arctg(2u),
where g1and g2 are signal functions given by (3.3)g1(pa,sign⁡(u))=βΔpa={(psup⁡-pa)βenchifu≥0,(pa-patm)βesvifu<0,g2(pb,sign⁡(u))=βΔpb={(psup⁡-pb)βenchifu<0,(pb-patm)βesvifu≥0,
where psup is the supply pressure, patm is the atmospheric pressure, and βench and βesv are the constant coefficients.Equations (3.2) are a fitting of a surface obtained experimentally [5, 7] in the test rig of Figure 1, considering that the piston is stopped; in that way, the volume is constant, and the speed of the piston is null. The mass flow rates at different pressures and valve input voltages were first estimated from the pressure versus time responses obtained for step inputs in valve voltage and a fixed piston position.The fitted mass flow rate in valve orifice,qma, is plotted versus input voltage and pressure difference in Figure 6.Fitted model of mass flow rate.
(a)(b)Rao and Bone [1] used a second-order bipolynomial equation to fit this function. In a similar way, Perondi [10] used a third-order polynomial one. Bobrow and McDonell [17] use a curve fit for the change in internal energy as a function of cylinder pressure which is quadratic in u. One of the greatest problems in these equations found in the literature is the difficulty in isolating the signal u, necessary when a control methodology that considers the nonlinear characteristics of the system is used. Equations to mass flow rate proposed by Ritter et al. [5] are innovations that have advantages as easiness of computational implementation and differentiation.


### 3.3. Pressure Dynamics

The cylinder used in this modeling is symmetric and without spindle. In mathematical modeling, the pressure changes in the chambers are obtained using energy conservation laws. Figure7 shows a schematic drawing of cylinder used.Figure 7
Cylinder’s schematic drawing.The relationship between the air mass flow rate and the pressure changes in the chambers is obtained using energy conservation laws. According to Perondi [10], the energy balance yields(3.4)qmaT-paCpdVadt=1γRddt(paVa),
where T is the air supply temperature, qma is the air mass flow rate into chamber A, pa is the absolute pressure in chamber A, Cp is the specific heat of the air at constant pressure, Cv is the specific heat of the air at constant volume, γ=Cp/Cv is the ratio between the specific heat values of the air, R is the universal gas constant, and V̇a=(dVa/dt) is the volumetric flow rate. Assuming that the mass flow rates are nonlinear functions of the servo valve control voltage (u) and of the cylinder pressures, that is, qma=qma(pa,u) and qmb=qmb(pb,u), the total volume of chamber A is given by(3.5)Va=Ay+Va0,
where A is the cylinder cross-sectional area, y is the piston position, and Va0 is the initial volume of air in the line and at the chamber A extremity, including the pipeline. The change rate for this volume is V̇a=Aẏ, where ẏ is the piston velocity.In this manner, calculating the derivative term in the right hand side of (3.4), and using Cp=(γR)/(γ-1), we can solve this equation to obtain(3.6)ṗa=-AγẏAy+Va0pa+RγTAy+Va0qma(pa,u).Similarly, for chamber B of the cylinder, we obtain(3.7)ṗb=AγẏVb0-Aypb-RγTVb0-Ayqmb(pb,u).Note that the pressure dynamics in chambers A and B,ṗa and ṗb, given by (3.6) and (3.7), are functions of piston position y, piston velocity ẏ, and mass flow rate qm. Figure 8 shows numerical simulation results presented in Endler [7] for constant mass flow rate input, where pa and pb as a function of time are illustrated. The pressure curves in the cylinder’s chambers for the control signal u=2 volts are observed, where the valve is opened, so that its “b” port is connected to the atmosphere.Figure 8
The numerical study of pressure dynamics behaviour in pneumatic cylinder’s chambers.More detailed study of (3.6) and (3.7) is presented in Bobrow and McDonell [17], Perondi [10] and Ritter et al. [5].


### 3.4. Friction Dynamics in Motion Equation

Applying Newton’s second law to the piston-load assembly results in(3.8)Mÿ+Fatr=Fp,
where M is the mass of the piston-load assembly, ÿ is the cylinder acceleration, Fatr is the friction force, and Fp is the pneumatic force related to the pressure difference between the two sides of the piston, that is given by A(pa-pb).In this section, the dynamic model to friction is based on the microscopic deformation of asperities in surface contact. It is possible to perceive an evolution in friction models that are based on the asperity microscopic deformations and depicted in recent papers.The Dahl model describes friction in the presliding movement phase, in a similar way to the rigid spring with damping behavior, but it has not included the Stribeck friction effect. The LuGre model, proposed by Canudas de Wit et al. [14], is an improved model that includes the Stribeck friction and describes many complex friction behaviors but is limited in the presliding movement phase, according to simulations results presented by Dupont et al. [15] and experimental tests carried out by Swevers et al. [18]. These authors also propose improvements in LuGre model through the inclusion of a model to hysteresis with nonlocal memory and sliding-force transition curves in presliding movement phase. This improved model is named Leuven model and used in friction modeling to a pneumatic servo positioning system by Nouri et al. [19]. Dupont et al. [15] also propose improvements in LuGre model through its interpretation as an elastoplastic friction model that is used in this paper.Figure9 represents the contact between surfaces through a lumped elastic asperity, considering a rigid body where the displacement y is decomposed into its elastic and plastic (inelastic) components z and w.Figure 9
Model of body subject to friction force showing elastic (z) and inelastic (w) displacement components.The friction force is described according to the LuGre friction model proposed by Canudas de Wit et al. [14]. In this model, the friction force is given by (3.9)Fatr=σ0ż+σ1z+σ2y,̇
where z is a friction internal state that describes the average elastic deflection of the contact surfaces during the stiction phases, σ0 is the stiffness coefficient of the microscopic deformations z during the presliding displacement, σ1 is a damping coefficient, σ2 represents the viscous friction, and ẏ is the velocity.The dynamicsż of the internal state z is modeled by the equation(3.10)dzdt=ẏ-α(z,ẏ)σ0gss(ẏ)|ẏ|z,
where gss(ẏ) is a positive function that describes the steady-state characteristics of the model for constant velocity motions and is given by(3.11)gss(ẏ)=Fc+(Fs-Fc)e-(ẏ/ẏs)2,
where Fc is the Coulomb friction force, Fs is the static friction force and ẏs is the Stribeck velocity. Figure 10 illustrate the behavior of the friction force as a function of velocity in steady state [8].Figure 10
Friction force characteristics combined in steady state.The functionα(z,ẏ) is presented according to Dupont et al. [15] and is used to represent the stiction. This function is defined by equations(3.12)α(z,ẏ)={0,if|z|≤zba0<12sen(πz-((zmax⁡(ẏ)+zba)/2)zmax⁡(ẏ)-zba)<1,ifzba<|z|<zmax⁡(ẏ)1,if|z|≥zmax⁡(ẏ)0,ifsign⁡(ẏ)≠sign⁡(z){0,if|z|≤zba0<12sen(πz-((zmax⁡(ẏ)+zba)/2)zmax⁡(ẏ)-zba)<1,ifzba<|z|<zmax⁡(ẏ)1,if|z|≥zmax⁡(ẏ)0,ifsign⁡(ẏ)≠sign⁡(z)},sign⁡(ẏ)=sign⁡(z),(3.13)0<zba<zmax⁡(ẏ)=gss(ẏ)σ0,para∀ẏ∈R,
where zba is a breakaway displacement, such a way that to z≤zba all movements in friction interface consist in elastic displacements only, and zmax is the maximum value of microscopic deformations and is velocity dependent.It is possible to observe that, withz represented by (3.12), when sliding movement is in steady state, ẏ is constant,  α(z,ẏ)=1, and ż=0. The z states values are given by equation(3.14)zss=ẏ|ẏ|gss(ẏ)σ0=sign⁡(ẏ)(Fc+(Fs-Fc)e-(ẏ/ẏs)2)σ0.Substituting (3.14) into (3.12), the friction force at steady state is obtained and is written as:(3.15)Fatrss=σ0zss+σ1⋅0+σ2ẏ=sign⁡(ẏ)(Fc+(Fs-Fc)e-(ẏ/ẏs)2)+σ2ẏ.A detailed study of (3.10), (3.11), and (3.15) is presented as follow. The properties of the friction model given by (3.9) and (3.10) will be explored. To capture the intuitive properties of the sliding body in Figure 9, the microscopic deformation z should be finite, that is, as important property, the state z is bounded.Proof.
Consider the following Lyapunov candidate function:(3.16)V=z22,
differentiating and combining with (3.10) and α(z,ẏ)=1= constant, it can be written as
(3.17)dVdt=-|ẏ|⋅|z|⋅(σ0⋅|z|gss(ẏ)-sign⁡(ẏ)⋅sign⁡(z)),
where dV/dt is negative if |z|>gss(ẏ)/σ0, since gss(ẏ) is strictly positive and bounded by Fs, see (3.11), then the set Ω={z:|z|≤Fs/σ0} is an invariant set for the solutions of (3.10).
Equation (3.15) is a nonlinear function that represents friction force as a function of velocity in steady state as was illustrated in Figure 10.
These dynamic properties of friction model presented are shown by Dupont et al. [15] and follow similar analysis carried out by Lyapunov method, as presented by Canudas de Wit et al. [14] and Canudas de Wit [20]. Among model main properties, it is cited that z state variable is limited and the model is dissipative, satisfies the stick and slip conditions, and represents adequately the presliding movement phase.
The applied force of Figure11(a) was chosen to challenge the stiction capability of the model; the force ramps up to cause breakaway and then returns to a level below that of Coulomb friction. Additionally, an oscillation is present and could be introduced by sensor noise or vibration. The response of friction model is seen in Figure 11(b). The friction dynamic model renders both presliding displacement and stiction.Applied force in pneumatic actuator and position response in both presliding displacement and stiction.
(a)(b)


## 3.1. Dead Zone Nonlinearity

This section presents the mathematical model for dead zone nonlinearity and its graphical representation. Dead zone is a static input-output relationship which for a range of input values gives no output. Figure4 shows a sectional view sketch of typical spool valve with main mechanical elements.Figure 4
Sectional view sketch of typical spool valve with main mechanical elements of the proportional valve with input dead zone.The mathematical model for dead zone in pneumatic servo valves presented in this section was obtained from Tao and Kokotović [16]. The dead zone analytical expression is given by the equation(3.1)uzm(t)={md(u(t)-zmd)ifu(t)≥zmd,0ifzme<u(t)<zmd,me(u(t)-zme)ifu(t)≤zme,
where u is the input value, uzm is the output value, zmd is the right limit of the dead zone, zme is the left limit of the dead zone, md is the right slope of the output, and me is the left slope of the output.Figure5 shows a graphical representation of the dead zone. In general, neither the break-points (zmd and zme) nor the slopes (md and me) are equal.Figure 5
Graphical representation of the dead zone.In current fluid power literature, dead zone in valves is expressed as a percentual of spool displacement. Valdiero et al. [6] present a new methodology for dead zone nonlinearity identification in proportional directional pneumatic valves. It is based on observing the dynamic behavior of the pressure in the valve gaps. The dead zone nonlinearity is a key factor that limits both static and dynamic performance in feedback control of fluid power systems. The usual method to cancel the harmful effects of dead zone is to add its fixed inverse function into the controller. This inverse is modeled by a set of parameters that need to be identified. The classic dead zone parameter identification uses expensive flow transducers and special test rig, while our proposed methodology needs only the pressure transducers shown in Figure 1. Experimental results are presented in Valdiero et al. [6] and illustrate the efficacy of this methodology that is cheaper and faster.


## 3.2. Mass Flow Rate

According to Rao and Bone [1], the mass flow rate model of the proportional valve is a key part of the system model. In this paper, we use an innovator model to mass flow rate equation qma and qmb developed by Endler [7], given by equations(3.2)qma(u,pa)=g1(pa,sign⁡(u))arctg(2u),qmb(u,pb)=g2(pb,sign⁡(u))arctg(2u),
where g1and g2 are signal functions given by (3.3)g1(pa,sign⁡(u))=βΔpa={(psup⁡-pa)βenchifu≥0,(pa-patm)βesvifu<0,g2(pb,sign⁡(u))=βΔpb={(psup⁡-pb)βenchifu<0,(pb-patm)βesvifu≥0,
where psup is the supply pressure, patm is the atmospheric pressure, and βench and βesv are the constant coefficients.Equations (3.2) are a fitting of a surface obtained experimentally [5, 7] in the test rig of Figure 1, considering that the piston is stopped; in that way, the volume is constant, and the speed of the piston is null. The mass flow rates at different pressures and valve input voltages were first estimated from the pressure versus time responses obtained for step inputs in valve voltage and a fixed piston position.The fitted mass flow rate in valve orifice,qma, is plotted versus input voltage and pressure difference in Figure 6.Fitted model of mass flow rate.
(a)(b)Rao and Bone [1] used a second-order bipolynomial equation to fit this function. In a similar way, Perondi [10] used a third-order polynomial one. Bobrow and McDonell [17] use a curve fit for the change in internal energy as a function of cylinder pressure which is quadratic in u. One of the greatest problems in these equations found in the literature is the difficulty in isolating the signal u, necessary when a control methodology that considers the nonlinear characteristics of the system is used. Equations to mass flow rate proposed by Ritter et al. [5] are innovations that have advantages as easiness of computational implementation and differentiation.


## 3.3. Pressure Dynamics

The cylinder used in this modeling is symmetric and without spindle. In mathematical modeling, the pressure changes in the chambers are obtained using energy conservation laws. Figure7 shows a schematic drawing of cylinder used.Figure 7
Cylinder’s schematic drawing.The relationship between the air mass flow rate and the pressure changes in the chambers is obtained using energy conservation laws. According to Perondi [10], the energy balance yields(3.4)qmaT-paCpdVadt=1γRddt(paVa),
where T is the air supply temperature, qma is the air mass flow rate into chamber A, pa is the absolute pressure in chamber A, Cp is the specific heat of the air at constant pressure, Cv is the specific heat of the air at constant volume, γ=Cp/Cv is the ratio between the specific heat values of the air, R is the universal gas constant, and V̇a=(dVa/dt) is the volumetric flow rate. Assuming that the mass flow rates are nonlinear functions of the servo valve control voltage (u) and of the cylinder pressures, that is, qma=qma(pa,u) and qmb=qmb(pb,u), the total volume of chamber A is given by(3.5)Va=Ay+Va0,
where A is the cylinder cross-sectional area, y is the piston position, and Va0 is the initial volume of air in the line and at the chamber A extremity, including the pipeline. The change rate for this volume is V̇a=Aẏ, where ẏ is the piston velocity.In this manner, calculating the derivative term in the right hand side of (3.4), and using Cp=(γR)/(γ-1), we can solve this equation to obtain(3.6)ṗa=-AγẏAy+Va0pa+RγTAy+Va0qma(pa,u).Similarly, for chamber B of the cylinder, we obtain(3.7)ṗb=AγẏVb0-Aypb-RγTVb0-Ayqmb(pb,u).Note that the pressure dynamics in chambers A and B,ṗa and ṗb, given by (3.6) and (3.7), are functions of piston position y, piston velocity ẏ, and mass flow rate qm. Figure 8 shows numerical simulation results presented in Endler [7] for constant mass flow rate input, where pa and pb as a function of time are illustrated. The pressure curves in the cylinder’s chambers for the control signal u=2 volts are observed, where the valve is opened, so that its “b” port is connected to the atmosphere.Figure 8
The numerical study of pressure dynamics behaviour in pneumatic cylinder’s chambers.More detailed study of (3.6) and (3.7) is presented in Bobrow and McDonell [17], Perondi [10] and Ritter et al. [5].


## 3.4. Friction Dynamics in Motion Equation

Applying Newton’s second law to the piston-load assembly results in(3.8)Mÿ+Fatr=Fp,
where M is the mass of the piston-load assembly, ÿ is the cylinder acceleration, Fatr is the friction force, and Fp is the pneumatic force related to the pressure difference between the two sides of the piston, that is given by A(pa-pb).In this section, the dynamic model to friction is based on the microscopic deformation of asperities in surface contact. It is possible to perceive an evolution in friction models that are based on the asperity microscopic deformations and depicted in recent papers.The Dahl model describes friction in the presliding movement phase, in a similar way to the rigid spring with damping behavior, but it has not included the Stribeck friction effect. The LuGre model, proposed by Canudas de Wit et al. [14], is an improved model that includes the Stribeck friction and describes many complex friction behaviors but is limited in the presliding movement phase, according to simulations results presented by Dupont et al. [15] and experimental tests carried out by Swevers et al. [18]. These authors also propose improvements in LuGre model through the inclusion of a model to hysteresis with nonlocal memory and sliding-force transition curves in presliding movement phase. This improved model is named Leuven model and used in friction modeling to a pneumatic servo positioning system by Nouri et al. [19]. Dupont et al. [15] also propose improvements in LuGre model through its interpretation as an elastoplastic friction model that is used in this paper.Figure9 represents the contact between surfaces through a lumped elastic asperity, considering a rigid body where the displacement y is decomposed into its elastic and plastic (inelastic) components z and w.Figure 9
Model of body subject to friction force showing elastic (z) and inelastic (w) displacement components.The friction force is described according to the LuGre friction model proposed by Canudas de Wit et al. [14]. In this model, the friction force is given by (3.9)Fatr=σ0ż+σ1z+σ2y,̇
where z is a friction internal state that describes the average elastic deflection of the contact surfaces during the stiction phases, σ0 is the stiffness coefficient of the microscopic deformations z during the presliding displacement, σ1 is a damping coefficient, σ2 represents the viscous friction, and ẏ is the velocity.The dynamicsż of the internal state z is modeled by the equation(3.10)dzdt=ẏ-α(z,ẏ)σ0gss(ẏ)|ẏ|z,
where gss(ẏ) is a positive function that describes the steady-state characteristics of the model for constant velocity motions and is given by(3.11)gss(ẏ)=Fc+(Fs-Fc)e-(ẏ/ẏs)2,
where Fc is the Coulomb friction force, Fs is the static friction force and ẏs is the Stribeck velocity. Figure 10 illustrate the behavior of the friction force as a function of velocity in steady state [8].Figure 10
Friction force characteristics combined in steady state.The functionα(z,ẏ) is presented according to Dupont et al. [15] and is used to represent the stiction. This function is defined by equations(3.12)α(z,ẏ)={0,if|z|≤zba0<12sen(πz-((zmax⁡(ẏ)+zba)/2)zmax⁡(ẏ)-zba)<1,ifzba<|z|<zmax⁡(ẏ)1,if|z|≥zmax⁡(ẏ)0,ifsign⁡(ẏ)≠sign⁡(z){0,if|z|≤zba0<12sen(πz-((zmax⁡(ẏ)+zba)/2)zmax⁡(ẏ)-zba)<1,ifzba<|z|<zmax⁡(ẏ)1,if|z|≥zmax⁡(ẏ)0,ifsign⁡(ẏ)≠sign⁡(z)},sign⁡(ẏ)=sign⁡(z),(3.13)0<zba<zmax⁡(ẏ)=gss(ẏ)σ0,para∀ẏ∈R,
where zba is a breakaway displacement, such a way that to z≤zba all movements in friction interface consist in elastic displacements only, and zmax is the maximum value of microscopic deformations and is velocity dependent.It is possible to observe that, withz represented by (3.12), when sliding movement is in steady state, ẏ is constant,  α(z,ẏ)=1, and ż=0. The z states values are given by equation(3.14)zss=ẏ|ẏ|gss(ẏ)σ0=sign⁡(ẏ)(Fc+(Fs-Fc)e-(ẏ/ẏs)2)σ0.Substituting (3.14) into (3.12), the friction force at steady state is obtained and is written as:(3.15)Fatrss=σ0zss+σ1⋅0+σ2ẏ=sign⁡(ẏ)(Fc+(Fs-Fc)e-(ẏ/ẏs)2)+σ2ẏ.A detailed study of (3.10), (3.11), and (3.15) is presented as follow. The properties of the friction model given by (3.9) and (3.10) will be explored. To capture the intuitive properties of the sliding body in Figure 9, the microscopic deformation z should be finite, that is, as important property, the state z is bounded.Proof.
Consider the following Lyapunov candidate function:(3.16)V=z22,
differentiating and combining with (3.10) and α(z,ẏ)=1= constant, it can be written as
(3.17)dVdt=-|ẏ|⋅|z|⋅(σ0⋅|z|gss(ẏ)-sign⁡(ẏ)⋅sign⁡(z)),
where dV/dt is negative if |z|>gss(ẏ)/σ0, since gss(ẏ) is strictly positive and bounded by Fs, see (3.11), then the set Ω={z:|z|≤Fs/σ0} is an invariant set for the solutions of (3.10).
Equation (3.15) is a nonlinear function that represents friction force as a function of velocity in steady state as was illustrated in Figure 10.
These dynamic properties of friction model presented are shown by Dupont et al. [15] and follow similar analysis carried out by Lyapunov method, as presented by Canudas de Wit et al. [14] and Canudas de Wit [20]. Among model main properties, it is cited that z state variable is limited and the model is dissipative, satisfies the stick and slip conditions, and represents adequately the presliding movement phase.
The applied force of Figure11(a) was chosen to challenge the stiction capability of the model; the force ramps up to cause breakaway and then returns to a level below that of Coulomb friction. Additionally, an oscillation is present and could be introduced by sensor noise or vibration. The response of friction model is seen in Figure 11(b). The friction dynamic model renders both presliding displacement and stiction.Applied force in pneumatic actuator and position response in both presliding displacement and stiction.
(a)(b)


## 4. Results

The most common and simple industrial application is a positioning task. By a positioning task, the objective of bringing the load position to a specified target in the actuator’s curse is meant. The proposed nonlinear mathematical model of a pneumatic servo position system was used in computer simulations of three cases of different cylinder size where the desired target position is 0.045 m.The pneumatic servo position system model dynamics is given by (3.1), (3.2), (3.6), (3.7), (3.8), (3.9), and (3.10). This model was implemented on the MatLab/Simulink software of which block diagram is shown in Figure 12 and using parameters presented in Tables 1 and 2. The classic proportional controller (P) was chosen because it is easy to implement and has only one parameter to adjust. Also, the results are easier to see with P controller.Table 1
The system parameters used in the numerical simulations.
Parameter valueDescriptionpsup=7×105 PaSupply pressurepatm=1×105 PaAtmospheric pressureR=287 Jkg/KGas constantT=293 KTemperature of supply airγ=1.4Specific heat ratio, dimensionlessM=10.125 kgMassβench=0.69501×10-8Constant coefficient to fill upβesv=0.898105×10-8Constant coefficient to deflateσ0=0.5×106 N/mStiffness coefficientσ1=400 Ns/mDamping coefficientσ2=46.532 Ns/mViscous friction coefficientFc=23.103 NCoulomb friction forceFs=24 NStatic friction forceyṡ=0.010 m/sStribeck velocityTable 2
The cylinder parameters used in the numerical simulations for different sizes.
Case (a)Case (b)Case (c)Descriptiond=0.025 m,d=0.060 m,d=0.060 mDiameterL=1 mL=1 mL=0.1746 mCurse lengthA=4.91×10-4 m2A=2.8×10-3 m2A=2.8×10-3 m2Piston areaVa0=2.5×10-4 m3Va0=1.4×10-3 m3Va0=2.5×10-4 m3Initial volume at chamber AVb0=2.5×10-4 m3Vb0=1.4×10-3 m3Vb0=2.5×10-4 m3Initial volume at chamber BFigure 12
Proportional feedback control structure block diagram with nonlinear mathematical model of pneumatic servo positioning system.The choice curse length for case (c) in Table2 is determined in such a way that it results in the same chamber volume of the case (a). It is a good idea because the chamber volume has a great influence in pressure dynamics given by (3.6) and (3.7). Results obtained from feedback control to positioning task simulations are depicted in Figure 13.Figure 13
Simulation results to positioning task: target positionyd=0.045 m, case (a) cylinder with small diameter, case (b) cylinder with large diameter, and case (c) cylinder with large diameter and small curse length.We are aiming at a good knowledge of the dynamic behavior for different pneumatic cylinders in order to see which features contribute to a better performance in a given engineering application. The results presented in case (a) outline the faster response with oscillating and overshoot in actuator position. Also, there are hunting problems that are oscillations caused by limit cycles around desired position. In many applications as robotics and aerospace engineering, the faster response is one of the requirements for the positioning task, and we can design pneumatic positioning systems with smaller cylinder diameter and increase the supply pressure obtaining necessary actuator force. To solve this overshoot problem, we can use an optimal control design for nonlinear systems as in the work of Rafikov et al. [12]. Also, the friction compensation is especially important, so that there are no hunting problems [4] and the actuator has an accurate response.Despite being slow, the results in case (b) are very good for some engineering applications as automatic welding, machining processes, surface finishing, and agricultural machinery, where application requirements do not permit overshoot, and task velocity is smaller. In this case, we can design pneumatic positioning systems with larger cylinder diameter that result in the damping increase, making the system slower. Besides, we can design a classical feedback control system depending on necessary accuracy in application. The case (c) presented dynamic behavior similar to case (b), and it shows that the chamber volume does not have significant influence in this positioning task.


## 5. Conclusion

This paper presented a full nonlinear mathematical model for pneumatic servo position system that can be used in numerical simulations to mechanical design and control system design of industrial applications. There was a bibliographical revision in recent international literature. However, these works do not address completely all the important nonlinearities in mathematical model. So, the main paper contribution was to present their nonlinearities and its complete mathematical modeling with some innovation and application results. The proposed systematic methodology is important to help researchers in the nonlinear modeling and precision control success. Future research will include an optimal nonlinear control strategy to overcome problems of the servo pneumatic system in agricultural machinery applications with high performance.



---

*Source: 472903-2011-05-23.xml*