Design, Modeling, Identification and Control of Mechatronic Systems

Oriol Gomis i Bellmunt Doctoral Thesis Barcelona, February 2007 DESIGN, MODELING, IDENTIFICATION AND CONTROL OF MECHATRONIC SYSTEMS Part I. Design Rules and Actuator Modeling for the Optimization of Mechatronic Systems Part II. Identification and Control of Piezoelectric Actuators Oriol Gomis i Bellmunt DESIGN, MODELING, IDENTIFICATION AND CONTROL OF MECHATRONIC SYSTEMS Departament d'Enginyeria Elèctrica Universitat Polit` ecnica de Catalunya Departament d’Enginyeria El` ectrica Doctoral Thesis Design, Modeling, Identification and Control of Mechatronic Systems Part I. Design Rules and Actuator Modeling for the Optimization of Mechatronic Systems Part II. Identification and Control of Piezoelectric Actuators Author: Oriol Gomis i Bellmunt Advisors: Samuel Galceran i Arellano Fay¸cal Ikhouane Barcelona, February 2007 ii S´on moltes les persones que m’han donat suport durant aquests anys i que han fet possible que ara estigui acabant d’escriure les ´ultimes l´ınies de la tesi doctoral. Moltes gr`acies a tots. Als directors de la tesi. Al Samuel Galceran pel seu suport incondicional, per comunicar-me la seva visi´o pr`actica de l’enginyeria i per fer f`acil el que podria haver estat dif´ıcil i al Fay¸cal Ikhouane per transmetre’m la seva metodologia cient´ıfica i la capacitat de tractar els problemes d’enginyeria des de la formalitzaci´o te`orica. A tots els companys d’Engitrol, en especial a l’Oscar i la Carmen, per introduir-me al m´on de l’enginyeria i a la realitat de la industria i pels anys que vaig passar amb ells. Den Kollegen und Freunden aus Braunschweig und vom DLR, mit denen ich immer interessante technische und nicht-technische Diskussionen hatte, die mir beim Deutsch halfen und mit denen ich eine tolle Zeit hatte. Flavio Campanile, f¨ur sein Rat und Beistand zum ersten Teil der Dissertation. Robert, Damiano, Bj¨orn, Nicole und Valerie, f¨ur ihre Freundschaft. A tots els companys del CITCEA i del Departament d’Enginyeria El`ectrica, pel suport que he rebut durant aquest temps i per l’intercanvi d’idees. Al Pere Castell, per la seva col·laboraci´o en els muntatges experimentals. Al Quim L´opez i al Daniel Montesinos per la seva valuosa ajuda en el disseny de les plaques i en la programaci´o del DSP. A l’Antoni Sudri`a per donar-me l’oportunitat d’integrar-me a l’equip del CITCEA. A la meva fam´ılia, per donar-me suport en tot moment. Als meus pares, a l’Anna, i als avis, en especial al padr´ı Crist`ofol, per transmetre’m la seva actitud positiva davant la vida. Als que heu aguantat les llargues explicacions quan em pregunt`aveu de qu`e anava aix`o de la tesi. A l’Andreu i el Xavi per tots aquests anys d’anades i vingudes. A l’Edu, el Josep, el Manuel i l’Oriol pels anys d’estudiants on potser no ens vam fer enginyers, per`o ens vam fer persones. Al Salva, el Joan i la Marta, per les llargues converses sense conclusions i per algunes conclusions sense converses. A la S´ılvia, per ser com ´es, per estar sempre al meu costat, per tot el que hem viscut i sobretot, per tot el que ens queda per viure. A la hist`eresis, la no-linealitat i la complexitat, perqu`e en aquest m´on on tot es vol simplificar i resumir en una frase, una de les poques coses clares que podem tenir ´es que les respostes s´on sempre dif´ıcils i complexes. La hist`eresis i les no-linealitats no s´on nom´es presents a l’enginyeria. Acknowledgements Supported by CICYT through grants DPI2005-08668-C03-03 and DPI2005-08668-C03-01. Part of this thesis is a result of a research training project which has been developed in the DLR (German Aerospace Center) in Braunschweig (Germany) and supported by a Marie Curie Fellowship of the European Community programme Smart Lightweight Structures And Transportation Application under the contract number HPMT-CT-2001-00298. vi Resum Les societats modernes plantegen nous reptes que demanden noves maneres de tractar els projectes d’enginyeria. Els enginyers han d’afrontar aquests reptes i desenvolupar solucions `optimes i eficients pels problemes cl`assics i nous. Els diferents aven¸cos produ¨ıts en la tecnologia hi poden ajudar, per`o una nova manera de tractar els problemes enginyerils ´es tamb´e necess`aria, no considerant ´unicament les diferent especialitats de l’enginyeria a¨ılladament. En aquest context, podem parlar de la creaci´o d’una nova filosofia de fer enginyeria: la Mecatr`onica. En els darrers anys han aparegut diferents definicions: A [2] la Mecatr`onica es defineix com l’aplicaci´o de decisions complexes a l’operaci´o de sistemes f´ısics. A [22;48] la Mecatr`onia ´es definida com la integraci´o o sinergia de diferents disciplines de l’enginyeria. Aquestes disciplines inclouen l’enginyeria mec`anica,l’enginyeria el`ectrica,l’enginyeria electr`onica,l’enginyeria de control,les comunicacions industrials il’enginyeria de software. A [50] es d´ona una definici´o m´es espec´ıfica: En general, la Mecatr`onica s´on soluciones de sistemes, que poden ser realitzades utilitzant components mec`anics, electr`onics, computacionals, materials, qu´ımics i de programari amb les seves corresponents disciplines enginyerils. L’objectiu d’aquestes solucions ´es incrementar la funcionalitat del sistema, l’intel·lig`encia i la fiabilitat, reduint els costos de producci´o. No obstant, tal com sost´e [48], la import`ancia del concepte no est`a ´unicament en la definici´o sin´o a la filosofia que hi ha al fons. ´ Es important de veure, que la Mecatr`onica no ´es nom´es la suma dels resultats de diferents disciplines, sin´o la filosofia enginyeril per xiv Thesis Outline The thesis has been divided into two parts. Part Iis entitled Design Rules and Actuator Modeling for the Optimization of Mechatronic Systems and includes chapters 1,2,3,4and 5.Part II is entitled Identification and Control of Piezoelectric Actuators and includes chapters 6, 7,8,9,10,11 and 12. The first part has been structured as follows. In chapter 1a brief introduction is presented, defining the orientation, motivation and objectives. In chapter 2the methodology is introduced, including the description of the different steps: design parameters, force-stroke and work-stroke curves, discussion of the limiting quantities involved in the output quantities expressions, maximum force for a given size, analysis of the scalability of the actuator, dimensional analysis study and comparison between the theoretical results and the industrial actuators quantities. In chapter 3the methodology is applied to electromagnetic actuators. Two linear electromagnetic actuators are considered: solenoids and moving coil actuators. In chapter 4the methodology is applied to hydraulic actuators. Finally, in chapter 5the conclusions are summarized. The second part has been structured as follows. In chapter 6a brief introduction is presented, defining the motivation and objectives. In chapter 7 an introduction to the piezoelectric effect is exposed. The relevant quantities and constants are introduced and the typical linear formulation for low and high frequency are stated. Piezoelectric actuators are introduced, along with their common applications, advantages and drawbacks. The Bouc-Wen hysteresis model is introduced in chapter 8. Thereafter, an identification technique to determine the parameters is proposed and its xv robustness against different classes of perturbations is discussed in chapter 9. Chapter 10 adds an adaptation to the previous model, in order to allow it to characterize better the behavior of piezoelectric actuators. The model is validated with a real actuator and the advantages over the previous chapter method are shown. In chapter 11 the models developed are employed to design a new controller. The controller takes into account not only the error but the output control effort which is tried to keep unchanged when a perturbation occurs. In chapter 12 the conclusions of the thesis part are summarized. xvi Contents I Design Rules and Actuator Modeling for the Optimization of Mechatronic Systems 1 1 Introduction 3 1.1 Design Rules and Actuator Modeling for the Optimization of MechatronicSystems....................... 3 1.2 Application of the methodology . . . . . . . . . . . . . . . . . 5 1.3 Objectives and Scope . . . . . . . . . . . . . . . . . . . . . . . 7 1.4 Outline............................... 7 2 Design Rules and Actuator Modeling for the Optimization of Mechatronic Systems 9 3 Application to electromagnetic actuators 15 3.1 Solenoidactuators ........................ 18 3.2 Moving coil actuators . . . . . . . . . . . . . . . . . . . . . . . 26 3.3 Industrial actuators . . . . . . . . . . . . . . . . . . . . . . . . 32 3.4 Conclusions ............................ 33 4 Application to hydraulic actuators 37 4.1 Step 1. Design parameters . . . . . . . . . . . . . . . . . . . . 38 4.2 Step 2. Force-stroke and work-stroke characteristic . . . . . . 40 4.3 Step 3. Limiting quantities. . . . . . . . . . . . . . . . . . . . 40 4.4 Step 4. Maximum force, stroke and work. . . . . . . . . . . . . 42 4.4.1 Forwardmotion...................... 42 4.4.2 Backward motion . . . . . . . . . . . . . . . . . . . . . 43 xvii CONTENTS 4.4.3 Considering forward and backward motion . . . . . . . 44 4.4.4 StrokeandWork ..................... 45 4.5 Step5.Scalability ........................ 48 4.6 Step 6. Dimensional Analysis . . . . . . . . . . . . . . . . . . 48 4.7 Step 7. Industrial actuators . . . . . . . . . . . . . . . . . . . 49 4.8 Conclusions ............................ 50 5 Conclusions 53 5.1 Contributions ........................... 53 5.2 Futurework............................ 54 II Identification and Control of Piezoelectric Actuators 55 6 Introduction 57 6.1 Modeling and validation of piezoelectric actuators . . . . . . . 58 6.2 Control of piezoelectric actuators considering the hysteresis . . 59 6.3 Objectives............................. 60 6.4 Outline............................... 61 7 Piezoelectricity 63 7.1 The piezoelectric effect . . . . . . . . . . . . . . . . . . . . . . 63 7.1.1 Abriefhistory ...................... 66 7.1.2 Deformation modes. . . . . . . . . . . . . . . . . . . . 68 7.2 Piezoelectric actuator simplified model . . . . . . . . . . . . . 69 7.2.1 Lowfrequency....................... 69 7.2.2 Highfrequency ...................... 70 7.2.3 Load............................ 72 7.2.3.1 Example..................... 72 7.3 Considerations .......................... 74 7.3.1 Non-linearities....................... 74 7.3.2 Temperature dependance . . . . . . . . . . . . . . . . . 75 7.3.3 Aging ........................... 76 xviii CONTENTS 7.3.4 Piezoelectric materials . . . . . . . . . . . . . . . . . . 76 7.4 Applications............................ 76 8 The Bouc-Wen model 81 8.1 Introduction............................ 81 8.2 The normalized Bouc-Wen model . . . . . . . . . . . . . . . . 82 8.2.1 Classification of the Bouc-Wen models . . . . . . . . . 82 8.2.2 The normalized Bouc-Wen model . . . . . . . . . . . . 84 9 Analysis and parameter identification of the Bouc-Wen model 87 9.1 Parameter identification for the Bouc-Wen model . . . . . . . 89 9.1.1 Classofinputs ...................... 89 9.1.2 Analytic description of the forced limit cycle for the Bouc-Wenmodel ..................... 90 9.1.3 Identification methodology . . . . . . . . . . . . . . . . 91 9.1.4 Robustness of the identification method . . . . . . . . 94 9.2 Numerical simulation example . . . . . . . . . . . . . . . . . . 98 9.3 Conclusion.............................101 10 Adaptation of the Bouc-Wen model for the modeling and validation of a piezoelectric actuator 103 10.1 Experimental observations . . . . . . . . . . . . . . . . . . . . 104 10.2 The modified model and the corresponding identification methodology................................105 10.2.1Modifiedmodel......................106 10.2.2 Non-hysteretic term parameter identification . . . . . . 106 10.2.3 Hysteretic term parameter identification . . . . . . . . 108 10.3 Piezoelectric actuator modeling . . . . . . . . . . . . . . . . . 109 10.3.1 Experimental setup . . . . . . . . . . . . . . . . . . . . 109 10.3.2 Identification procedure . . . . . . . . . . . . . . . . . 111 10.3.3 Model validation . . . . . . . . . . . . . . . . . . . . . 113 10.4Conclusion.............................115 xix CONTENTS 11 Control of a piezoelectric actuator considering the hysteresis121 11.1 Background results. PID control of a Bouc-Wen hysteresis . . 122 11.2 Experimental Platform . . . . . . . . . . . . . . . . . . . . . 124 11.2.1 Experimental Layout . . . . . . . . . . . . . . . . . . . 124 11.2.2 System modeling . . . . . . . . . . . . . . . . . . . . . 127 11.2.3 Control objective . . . . . . . . . . . . . . . . . . . . . 127 11.3 Parameter identification . . . . . . . . . . . . . . . . . . . . . 128 11.4Controllaws............................131 11.4.1PIDControl........................131 11.4.2 PID plus a sinusoidal component . . . . . . . . . . . . 134 11.4.3 PID plus a sinusoidal component with a time varying amplitude .........................137 11.5 Experimental Results . . . . . . . . . . . . . . . . . . . . . . . 139 11.5.1PIDControl........................139 11.5.2 PID plus a sinusoidal component . . . . . . . . . . . . 141 11.5.3 PID plus a sinusoidal component with a time varying amplitude .........................141 11.6Conclusion.............................141 12 Conclusions 145 12.1Contributions ...........................145 12.2Futurework............................146 References 153 A Outline of the proof of Theorem 3 155 A.1 Determination of the parameter κx...............155 A.2 Existence and unicity of the zero of the function θ◦(¯x).....156 A.3 Determination of the parameter n................157 A.4 Determination of the parameter κw...............158 A.5 Determination of the parameter ρ................159 A.6 Determination of the parameter σ................159 xx CONTENTS B Publications 161 B.1 Journalpapers ..........................161 B.1.1 Published .........................161 B.1.2 Submitted.........................161 B.2 Conferencepapers.........................162 B.2.1 Published .........................162 B.2.2 Accepted..........................162 xxi CONTENTS xxii List of Tables 3.1 Solenoid actuator dimensional analysis quantities . . . . . . . 25 3.2 Moving coil actuator dimensional analysis quantities . . . . . . 31 4.1 Hydraulic actuator dimensional force analysis quantities . . . . 48 4.2 Hydraulic actuator dimensional work analysis quantities . . . 48 7.1 Piezoelectric material relevant parameters. . . . . . . . . . . . 76 7.2 Piezoelectric material properties [52]............... 77 7.3 Main applications of piezoelectric devices. . . . . . . . . . . . 79 8.1 Classification of the BIBO, passive and thermodynamically consistent Bouc-Wen models . . . . . . . . . . . . . . . . . . . 83 8.2 Classification of the BIBO, passive and thermodynamically stable normalized Bouc-Wen models . . . . . . . . . . . . . . . 85 10.1 Design parameters. . . . . . . . . . . . . . . . . . . . . . . . 109 10.2 Parameter expressions . . . . . . . . . . . . . . . . . . . . . . 109 10.3 The gicoefficients.........................113 10.4 Coefficients κi...........................113 10.5 Bouc-Wen model parameters . . . . . . . . . . . . . . . . . . . 114 11.1 Identified parameters. . . . . . . . . . . . . . . . . . . . . . . . 128 xxiii Chapter 1 Introduction An actuator can be defined as an energy converter which transforms energy from an external source into mechanical energy in a controllable way. The actuator input quantities depend on the type of energy used. For electromagnetical, piezoelectric and magnetostrictive actuators the input quantities can be the current, the charge or the voltage; for fluid power actuators the fluid pressure or the flow; for shape memory alloys and thermal expansion actuators the temperature. The relevant output quantities to be considered in the optimization are the force, the work and the stroke. The input quantities are provided by a control system which lead output quantities to the referenced values. Such quantities are ruled by the mechanical load system or structure, which defines the relationship between the force and the stroke. The integration of actuators and loads in a mechatronic or adaptronic1system allows the conception of a unique system which is to be analyzed. 1.1 Design Rules and Actuator Modeling for the Optimization of Mechatronic Systems The increasing quantity of different novel actuator technologies being used in different industrial applications along with the need for light and volume 1Adaptronics is a term referred to the analysis, design and integration of smart structures and systems. 3 1. Introduction reduced systems are boosting the necessity of general analysis using uniform criteria. Regarding the comparison of different actuators, in [1] an actuator selection criterium is presented in order to develop a software to choose the most suitable actuator for different applications. To undertake such a task performance indices and material property charts are provided. The actuators as energy converters are analyzed and compared in [37], focusing on robot applications. A comparison of the performance of different actuators regarding stress, strain, energy and precision is introduced in [23], showing different tables and graphics to compare the performance of the studied actuators. In [12] the performance of solid-state actuators available in the market is compared and studied. The environmental impact of the mechanical design is presented in [20], and material property charts including this new criteria are presented. A new selection and classification criterium is introduced in [60], including a comparison between existing actuators in the market as well as their stress, strain, power densities, and resolution. Although the previous references deal with the comparison of different actuators using different criteria, to the best of our knowledge no methodology to address the modeling of each actuator allowing comparisons of different classes of actuators has been found. Therefore, the present work proposes a novel methodology which might be applied to any class of actuator in order to optimize mechatronic and adaptronic systems. The motivation to develop a methodology stems from the need for light and volume reduced structures and systems, which are to be integrated in the design procedure as early as possible. Hence, the geometric relationships, aspect ratios and material properties that maximize the actuator output quantities with a certain limited volume or weight, along with their scalability for the integration in structures are studied. A validation of the results is done by performing dimensional analysis of the expressions obtained and comparing numerical results with industrial actuator data. 4 1.2 Application of the methodology 1.2 Application of the methodology Once the methodology has been presented, it is applied to electromagnetic and hydraulic actuators. Electromagnetic actuators are commonly used in many engineering fields. They have good force and work densities, although not as high as hydraulic actuators. They are easily controllable and the power source providing the energy can be placed as far away as necessary. Their use must be avoided when their environment must be free of electromagnetic fields or interferences. However, many technologies to deal with such effects are being developed. The electrical circuit provides the current to the coils. This current flows through wires and produces heat due to the Joule effect. Different materials can be employed for the wires but usually copper, silver or aluminium are used because they present the lowest resistivities. The magnetic circuit provides the flux and the force. Different materials can be used in the magnetic circuit depending on the magnetic permittivity. The electromechanical actuators have an electrical and a magnetic circuit. Such circuits are built together, and therefore, the heat generated in the coils by Joule effect must flow through part of the magnetic circuit. The heat transfer circuit includes all the components of the actuator and depends on the geometry of each of them. Although different materials can be used in both the electric and magnetic circuit, it is common to talk about copper for the electric circuit and iron for the magnetic circuit. As far as fluid power actuators are concerned, they use the fluid power to provide mechanical work; the difference between the pressures Pin two different chambers results in a relative pressure which produces a force Fin a given surface Swhich yields F=PS. The pressure is the input quantity, performing the same function as the current in electromechanical actuators. The fluid actuators employed in the industry are mainly divided by the state of the fluid employed: hydraulic actuators employ an incompressible liquid (usually oil), while pneumatic actuators employ a compressible gas (air). Hydraulic actuators are commonly used in many engineering fields. They show the following advantages: very good force and work densities (more than any other actuator), strokes as long as necessary (if enough fluid is supplied), 5 1. Introduction easily controllable and the fact that the power source providing the energy can be placed far away from the actuator (but not as far as with the electromagnetic actuators). Their main disadvantages are the safety problems generated by the high pressures needed (the same fact that provides the advantages), the leakage flow (that can become an important problem for actuator performance, safety conditions and environmental issues) and the hardly inflammability of the oil employed. Pneumatic actuators are used in many engineering fields, as well. They present good force and work densities, even though not as high as the hydraulic actuators, they can perform strokes as long as needed like their hydraulic counterparts, they are easily controllable and the power source providing the energy can be placed far away from the actuator. However, they cannot work with pressures as high as the hydraulic actuators because of the problems derived from the high compressibility of the gases. This same fact makes the hydraulic actuators faster in response and stiffer against external load disturbances. The efficiency of the hydraulic systems is also higher. It is caused by the losses of energy due to the heat transfer (in the air cooling), higher leakage and worse lubrication which occurs in the pneumatic systems. Nevertheless, the pneumatic systems can work at higher environment temperatures. In the present work, the hydraulic actuators are studied. However, some of the results obtained also apply for their pneumatic counterparts. Although several studies [1;12;23;37;60] in the literature provide rules and charts for selecting the most optimal actuator class for different applications and others [21] delve into the study of different actuator classes, to the best of our knowledge no detailed analysis of linear electromagnetic or hydraulic actuators following a general procedure and oriented towards improving the actuator design has been found. The present thesis introduces a new methodology to analyze linear electromagnetical and hydraulic actuators by modeling their maximum output mechanical quantities (force, work and stroke) as functions of the geometry and material properties and discusses the scalability (in the sense of producing the same stress and strain distribution for different sizes). 6 1.3 Objectives and Scope The motivation of this thesis part is to provide the detailed analysis of different actuators using a general procedure, rather than introducing a general analysis (see [23]). The static behavior1of different classes of actuators is considered. Such actuators include linear hydraulic and electromagnetic actuators. In this second group, the actuators have been chosen as examples of electromagnetic actuators with (moving coil) and without (solenoid) permanent magnets. 1.3 Objectives and Scope The objectives of the present thesis part may be summarized as: 1. Design of a methodology to deal with the modeling and optimization of industrial actuators. The methodology is the crucial step for the study of different classes of actuators. The methodology includes design optimization, scalability analysis and validation with real actuators and dimensional analysis. 2. Application of the methodology. The methodology has to be applied to common industrial actuators: hydraulic and electromagnetic. These two classes of actuators are the most employed in the industry. As far as the scope is concerned, the optimization is performed for linear hydraulic and electromagnetic actuators considering static behavior. Neither the dynamics nor the non-linear motion are considered in this work. 1.4 Outline The thesis part has been structured as follows. In chapter 2the methodology is introduced, including the description of the different steps: design parameters, force-stroke and work-stroke curves, discussion of the limiting quantities involved in the output quantities expressions, maximum force for 1The static behavior analysis assumes very slow operation, and hence it does not take into account the effect of the frequency in the analyzed system. 7 1. Introduction a given size, analysis of the scalability of the actuator, dimensional analysis study and comparison between the theoretical results and the industrial actuators quantities. In chapter 3the methodology is applied to electromagnetic actuators. Two linear electromagnetic actuators are considered: solenoids and moving coil actuators. In chapter 4the methodology is applied to hydraulic actuators. Finally, in chapter 5the conclusions are summarized. 8 Chapter 2 Design Rules and Actuator Modeling for the Optimization of Mechatronic Systems As previously stated, the purpose of this work is to develop design rules and models for actuator optimization. This section explains the general procedure introducing all the concepts which are going to provide such rules. The main steps are: 1. Design parameters. Study of the geometry and materials of the actuators. 2. Force-stroke and work-stroke curves. Analysis of the force, stroke and work production. 3. Discussion of the limiting quantities involved in the output quantities expressions. 4. Maximum force for a given size. Study of the limit force, stroke and work. 5. Analysis of the scalability of the actuator. 6. Dimensional analysis study of the relevant quantities. 9 2. Design Rules and Actuator Modeling for the Optimization of Mechatronic Systems 7. Comparison between the theoretical results and the industrial actuators quantities. The first step introduces the design parameters in the construction of the actuator. A detailed schematic drawing is presented showing the geometric properties and the materials used. In order to obtain a clear design parametrization geometrical factors,aspect ratios and filling factors are presented. The geometrical factors define the ratio between any geometrical dimension and a reference geometrical dimension in the same axis. A nondimensional factor kiis obtained for each length lias a quotient of this length and the reference length lin its axis as: ki=li l→li=kil(2.1) Using these factors, all the lengths in the same axis are related to one single length, simplifying the analysis of the size dependence of different quantities. The number nof independent reference lengths depend on the degrees of symmetry of the actuator. An actuator with cylindrical shape presents two different reference lengths (n= 2, since a cross-section diameter and a length define a cylinder), a spherical actuator would be defined with one reference length (n= 1, only a diameter defines a sphere). The relationships between different reference lengths is obtained using aspect ratios as: η=r l(2.2) If nindependent reference dimensions are necessary, n−1 aspect ratios are to be defined. The combination of the previous two concepts implies that all the geometric dimensions are expressed as a function of one single reference length, which is associated to the size of the actuator and allows the independent study of the performance of an actuator with a limited size and the actuator performance when the size is changed. The filling factor provides the portion of usable cross-section surface when electric wires are concerned. Due to the shape of the wires and the necessary 10 electrical isolation, the entire cross-section designed for the copper wires is not employed. The filling factor yields: kff =Scopper Stotal (2.3) where Scopper is the cross-section of the wires and Stotal the overall crosssection of the coil. A general expression of the output quantities as a function of all the input quantities involved is developed in the second step. These expressions are taken from the general physics laws ruling the actuators concerned. Each type of actuator behaves in a different way and its expressions are presented describing all the assumptions done. The force-stroke curves and work-stroke curves are obtained, establishing the characterization plot of the actuator. The inputs (currents, voltages, pressures, etc.) capable of changing these curves are presented, explaining why and how they can influence the actuator performance. The different working points depending on the load are discussed and graphically shown. It is important to note that the design parameters cannot be considered inputs and their influence is discussed in the following steps. The third step focuses on the quantities involved in the expressions obtained in the second step. The output quantities developed by an actuator can be controlled by modifying the input quantities (electrical voltages and currents, fluid pressures and flows, etc.). Some physical limits (maximum allowed temperature, mechanical resistance etc.) do not allow the actuator output quantities to be increased indefinitely. Since the purpose of the present work is to separately deal with the maximum force, stroke and work available in a given size and the performance scalability, only geometric quantities (reference lengths), relationships (geometrical factors, aspect ratios and filling factors), material properties (magnetic permeability, resistivity, resistivity temperature coefficient, conductivity, etc.), universal physics constants (µ0,σ, etc.) and physical thresholds (maximum temperature, stress, etc.) are to be used. Therefore, all the other quantities (currents, magnetic fluxes, pressures, etc.) must be expressed as functions of the mentioned quantities. 11 3. Application to electromagnetic actuators tion (since the air is surrounding de coils), the maximum current can be expressed as: imax =sAw∆Tπ(r2 out −r2 in)hc δ0(1 + γ∆T)lw (3.7) The conduction coefficient λis a material property, but the convection coefficient hcdepends on the non-dimensional Nusselt number which is expressed in [3] as: Nu−L=hcL λair (3.8) The Nusselt number can be written as a function of Reynolds, Prandtl and Grashof numbers, in [33] it is presented as: Nu=CRm ePn rGp r(3.9) where Reis the Reynolds number (ρvL/η) which shows the relationship between the inertial forces and the viscous forces in the dynamics of a fluid, Pris the Prandtl number (ηc/λ) which characterizes the regime of convection, Gris the Grashof number(βg∆TL3/ν2) analog to the Reynolds number when natural convection is concerned and C,m,nand pcan take different values in forced convection (C < 1, m < 1, n = 1/3, p = 0) and natural convection (C < 1, m = 0, n < 1/3, p < 1/3). The Nusselt number can be in all the cases expressed as: Nu−r=KNurα(3.10) where KNu and αmust be discussed in each case. If it is not otherwise stated the values used in the numerical calculations done in this work are λair = 0.0257 W/Km,λiron = 80 W/Km,µ0= 1.25664· 10−6Tm/At, ∆T= 50 K,µr= 200, µc=µm= 1, Hc= 0.5·106A/m, ρ0= 1.68 ·10−8Ωmand γ= 0.0068 Ωm/K. 3.1 Solenoid actuators Step 1.The solenoid actuators provide motion exciting a magnetic field where a plunger (movable part) tries to minimize the reluctance (i.e. the 18 3.1 Solenoid actuators air gap) moving to the less reluctance position. The geometry is shown in Fig. 3.1. The non-dimensional constant kψrefer to ψgeometry dimensions of Fig. 3.1. l2 l l1rr3r1 r2 hcu x copper wires iron pipe iron plates Figure 3.1: Solenoid actuator sketch Step 2.The magnetic flux flowing inside a solenoid can be derived from the reluctance expression. It can be written as: Φ = Fmm <=Ni x µ0S+l2+leq−x µrµ0S =Niµrµ0S l2+leq +x(µr−1) (3.11) where <is the reluctance expressed as a function of the magnetic properties of the iron µr, the length l2, the cross-section of the plunger S=πr2 1and the length leq, which is the plunger length with a reluctance equivalent to the reluctance of the plates and the pipe. Fmm is the magneto-motive force, equal to the number of turn Ntimes the current i. The number of turns can 19 3. Application to electromagnetic actuators be expressed as a function of the actuator dimensions as: N=hcul2kff Aw (3.12) where hcu is the thickness of copper, l2the coil’s length, kff the filling factor described in (2.3) and Awthe cross-section of a single wire. The solenoid force is produced for the change of reluctance due to the change of the air gap distance. Its expression can be derived from the energy stored in a solenoid Wm=Ridλ =RNidΦ. It yields: F=dWm dx =SN2i2µ2 rµ0 2 (l2+leq + (µr−1) x)2(3.13) The energy can be obtained integrating the force (3.13) between a given displacement xand 0 as: W=Zx 0 Fdx =SN2i2µ2 rµ0x 2 (l2+leq + (µr−1) x) (l2+leq)(3.14) It can be seen that Wis the total energy which the actuator stores in each position. This energy is transformed in work against a load and kinetic energy Wk= (1/2)mv2, since this work focuses on the static behavior of the studied actuators, a quasi-static movement is considered. Therefore, if it is not otherwise stated all the energy is assumed to be transformed into work. The force-strokes curves can be seen in Fig. 3.2. The force-stroke curves are presented when the input quantity (electrical current) is changed for different loads (one elastic load, equivalent to a structure, and one constant load). It can be seen how the operating points are changing depending on the input quantity and the load, when the current is increased working against an elastic load the working point is moving from E0 to E5. From the initial working point E0 the load can be moved to the other points depending on the input current. If a current 3i0is applied the plunger moves from E0 to E3 and remains there. The work against a constant load presents more difficulties. When the current is not large enough the actuator cannot begin to move and remains blocked at the initial posi- 20 3.1 Solenoid actuators 0 1 2 3 x 10−4 0 5000 10000 15000 Stroke Force o o o o o i=5i0 i=4i0 i=3i0 i=2i0 i=i0 Elastic Load Constant Loa d E0 E1 E2 E3 E4 S3 S4 S5 E5 o Figure 3.2: Force-Displacement curves for elastic and constant loads tion, for this reason S1 and S2 do not appear in the curve and a minimum current 2.5i0must be established to begin the traction. The force-stroke curves present hyperbolic behavior, with a horizontal asymptote tending to 0 as the distance increases and a vertical asymptote located in the negative stroke segment (physically unreachable). Step 3.The general expressions obtained for electromagnetic actuators (3.13), (3.12) and (3.8) apply for solenoid actuators. Step 4.Replacing the maximum current (3.6) in (3.13), using the number of turns given in (3.12), the convection coefficient of (3.8) and the geometrical expressions of (2.1) the maximum force (obtained when x= 0) can be expressed as: Fmax Sact =λairµ0∆Tµ2 r δ0(1 + γ∆T) k2 Lk2 r1kr3−kr1 kr2kff 2kλλair λiron +4 NuD(3.15) 21 3. Application to electromagnetic actuators with kλ=log(1/kr3) and kL=l2/(l2+leq). The equivalent length ratio can be expressed as: leq l=k2 r1(1 −2kl1) 1−k2 r3 +k2 r1η2log 1 kr1 kl1 (3.16) It can be noted that the maximum force divided into the cross-section of the actuator is expressed as a function of material constants, physical thresholds and geometrical relationships. A design factor depending on the design parameters can be defined from (3.15) as: qf=k2 Lk2 r1kr3−kr1 kr2kff (3.17) Substituting all the terms in the last expression it can be written as: qf=k2 l22k2 r1kff kr3−kr1 kr3+kr1 kl2+k2 r1kl2 1−k2 r3+k2 r1η2log 1 kr1 (1−kl2)/22(3.18) The expression (3.18) has been analyzed numerically. The best design parametrization has been found for values kr1= 0.34, kr3= 0.76, kl2= 0.50 and η < 0.1. The optimized found design factor is qf= 0.0402. In Fig. 3.3 the design factor qfdepending on the ratios kr1and kr3is plotted. The aspect ratio η, the ratio rl1and the filling factor are kept constant to allow a three-dimensional plot. The filling factor kff is typically around 0.75 and can be considered independent of the other design parameters. Regarding the aspect ratio, in (3.18) it is shown that a large ηimplies a low design factor qf, but its importance depends on the other terms on the denominator. It has been seen that below aspect ratios of 0.1 the improvement of the design factor is insignificant. In Fig. 3.3 it can be seen that kr1values between 0.3 and 0.4 provide the best performance for kr3values between 0.7 and 0.8. In Fig. 3.4 the design factor behavior depending on kl2and kr3with a constant kr1= 0.34 is presented. The plot shows that high design factors are obtained for kr3values between 0.7 and 0.8 and in a wide range (0.1−0.9) of kl2. 22 3.1 Solenoid actuators 0 0.5 1 0 0.2 0.4 0.6 0.8 1 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 kr3 Design factor kr1 Figure 3.3: Solenoid design factor depending on kr1and kr3with kl2= 0.5 and η= 0.1 The maximum displacement is l2and is proportional to the length of the actuator lsince l2=kl2l. The maximum volumetric work is achieved when the whole displacement is done. It can be obtained integrating the force. The maximum work expression found yields: Wmax Vact =qfλairµ0∆Tµ2 r kL (1/kL+µr−1) δ0(1 + γ∆T)2kλλair λiron +4 NuD(3.19) The discussion undertaken for the force optimum design parametrization applies for the work, too. Step 5.If the leq/l ratio is kept constant in (3.15) and (3.19) the scalability will depend only on the Nusselt number. If the Nusselt number is assumed to be constant (α= 0), the force will be independent of the actuator length and it will be proportional to the cross-section. The work will be proportional to both the length and the cross-section, and therefore to the volume so that a constant volumetric energy will be shown. Such an assump- 23 3. Application to electromagnetic actuators 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 kl2 Design factor kr3 Figure 3.4: Solenoid design factor depending on kl2and kr3with kr1= 0.34 and η= 0.1 tion cannot generally be done when studying the different convection cases. With a positive value of α(which is the behavior observed) the actuator force and work are not scalable anymore and their maximum performance values are increased when the size is increased. In such a case the maximum force per cross-section can be expressed as: Fmax Sact =ka rα kb+rα(3.20) It can be noted that for high values of rthe maximum force tends to be scalable since lim r→∞ karα/(kb+rα) = ka. It can also be observed in Fig. 3.5. For tiny actuators the force and work are strongly unscalable and the performance becomes worse. It may explain that these actuators are not used when a small actuator is required. Step 6.The force provided by these actuators can be analyzed with dimensional analysis using the Buckingham Pi Theorem [5], the quantities involved are shown in the Table 3.1. The FLTIθ (force - length - time - 24 3.1 Solenoid actuators 10−10 10−5 1001051010 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1Force scalability Radius Force / Cross−Section F/S α=0 α=0.25 α=0.5 α=0.75 α=1 Figure 3.5: Solenoid force scalability for different αcoefficients in the Nusselt number expression current - temperature) system is used. Table 3.1: Solenoid actuator dimensional analysis quantities FForce [F] xPosition [L] µPermeability [FI−2] hcConvection Coef. [FL−1T−1θ−1] λConduction Coef. [FT−1θ−1] δResistivity [FL2I−2T−1] TTemperature [θ] The results show that the maximum force can be expressed as: F x2=KµTλ δΦ (Nu) (3.21) where it can be observed that it matches perfectly with (3.15) since δ= δ0(1 + γ∆T) and µrand all the geometrical constants are adimensional. The results obtained for the work match with the expression (3.19), the maximum 25 3. Application to electromagnetic actuators work given by dimensional analysis can be expressed as: W x3=KµTλ δΦ (Nu) (3.22) 3.2 Moving coil actuators Step 1.The geometric dimensions are shown in Fig. 3.6. For the sake of simplicity and without loss of generality it has been assumed that l1=l2=l3. Step 2.Moving coil actuators use the force produced by the interaction of perpendicular magnetic field and electrical current, described in the Lorenz force law. It yields: F=Blwi(3.23) where Bis the field density provided by the permanent magnet, lwis the length of the wire and ithe current flowing in the wire. The present work assumes that the moving coil stroke is limited to the region where the flux is flowing, so that the force-stroke curve presents a constant force depending linearly on the current applied to the coil. The work is obtained by the integration of a constant function. Without the assumption of limited stroke, the length of wire lwbeing crossed by magnetic flux decreases as the coil is moving outside the flux region, whereas the flux density and the current are kept constant since the copper permittivity can be considered equal to the air permittivity. The flux density Bin the coil can be derived from the reluctance expression. In this case the magneto-motive force is provided by the permanent magnet Fmm =Hcl, where Hcis the coercitive field (a magnet constant). The reluctance can be calculated as the series association of all the reluctances sketched in Fig. 3.7. The total reluctance of the magnetic circuit can be expressed as: <= log k 1 kl1 −µr µcukl4 r2k µr−1 kl1 r4 k 2 kl1 r1k µr kl1 −µr µcukl4 r3!+2(1−2kl1) (1−k2 r5)η2+2 µm µrk2 r1η2 µrµ02πl (3.24) 26 3.2 Moving coil actuators L r1 r l magnetic circuit permanent magnet r2 r5 r3r4 d l2 l1 l3 l4 moving coil Figure 3.6: Moving coil actuator sketch 27 3. Application to electromagnetic actuators 10−8 10−7 10−6 10−5 10−4 10−3 10−2 10−1 100 10−1 100 101 102 103 104 105 106 107 Area [m2] Force [N] 1e+009 Pa 1e+007 Pa 1e+005 Pa Solenoid Ledex Solenoid Densitron Solenoid NSF Moving coil Beikimco Moving coil Maccon Maximum Solenoid Maximum Moving coil Regression Solenoid Regression Moving coil Figure 3.10: Industrial electromagnetic actuator force-area comparison 10−10 10−8 10−6 10−4 10−2 100 10−4 10−3 10−2 10−1 100 101 102 103 104 105 106 Volum [m3] Work [J] 1e+008 J/m3 1e+006 J/m3 1e+004 J/m3 Solenoid Ledex Solenoid Densitron Solenoid NSF Moving coil Beikimco Moving coil Maccon Maximum Solenoid Maximum Moving coil Regression Solenoid Regression Moving coil Figure 3.11: Industrial electromagnetic actuator work-volume comparison 34 3.4 Conclusions overcome by using certain design techniques including some high reluctance parts which are beyond the scope of this work. The limit output quantities have been found for certain aspect ratios and geometric relationships. The results have been compared with the performance of industrial actuators and it has been noted that the industrial actuators behave as expected. Regressions linking the force and work with the cross-section and volume have been carried out, resulting in similar performances as a function of the size as the theoretically developed models. The results can be used in any application with volume and weight constraints. For given volume or weight constraints, the presented expressions can show whether the studied electromagnetic actuators match the requirements and what materials and shapes are needed. 35 3. Application to electromagnetic actuators 36 Chapter 4 Application to hydraulic actuators Fluid power actuators use the fluid power to provide mechanical work; the difference between the pressures Pin two different chambers results in a relative pressure which produces a force Fin a given surface Swhich yields F=PS. The pressure is the input quantity, performing the same function as the current in electromechanical actuators. The fluid actuators employed in the industry are mainly divided by the state of the fluid employed: hydraulic actuators employ an incompressible liquid (usually oil), while pneumatic actuators employ a compressible gas (air). Hydraulic actuators are commonly used in many engineering fields. They show the following advantages: very good force and work densities (more than any other actuator), strokes as long as necessary (if enough fluid is supplied), easily controllable and the fact that the power source providing the energy can be placed far away from the actuator (but not as far as with the electromagnetic actuators). Their main disadvantages are the safety problems generated by the high pressures needed (the same fact that provides the advantages), the leakage flow (that can become an important problem for actuator performance, safety conditions and environmental issues) and the hardly inflammability of the oil employed. Pneumatic actuators are also used in many engineering fields. They present good force and work densities, even though not as high as the hydraulic actuators, they can perform strokes 37 4. Application to hydraulic actuators as long as needed like their hydraulic counterparts, they are easily controllable and the power source providing the energy can be placed far away from the actuator. However, they cannot work with pressures as high as the hydraulic actuators because of the problems derived from the high compressibility of the gases. This same fact makes the hydraulic actuators faster in response and stiffer against external load disturbances. The efficiency of the hydraulic systems is also higher. It is caused by the losses of energy due to the heat transfer (in the air cooling), higher leakage and worse lubrication which occurs in the pneumatic systems. Nevertheless, the pneumatic systems can work at higher environment temperatures. In the present work, the hydraulic actuators are studied. However, some of the results obtained also apply for their pneumatic counterparts. In this work an ideal power supply with no losses will be considered, it implies that the load will not change the supplied pressure and it can be assumed with no loss of generality if the power of the power supply is larger than the nominal power consumed by the actuator. The methodology described in chapter 2is applied below. 4.1 Step 1. Design parameters In Figure 4.1 a hydraulic actuator is sketched. It can be seen that for x= 0 both orifices are completely closed, when x > 0 follows P1> P2since Ps> Pr, and the plunger moves forward, when x < 0 follows P2< P1 and it moves backward. The sections can be written as A1=πD2 1/4 and A2=π/4 (D2 1−D2 2) where D1is the diameter of the cylinder and D2is the diameter of the rod which guides the plunger. The geometry of hydraulic cylindrical actuators is shown in Figure 4.2. The same geometry would be valid for pneumatic actuators, with the only difference of the fluid used and the corresponding limitations. 38 4.1 Step 1. Design parameters Figure 4.1: Hydraulic Actuator Figure 4.2: Geometry of a hydraulic actuator 39 4. Application to hydraulic actuators 4.2 Step 2. Force-stroke and work-stroke characteristic The cylinder force can be expressed as F=P1A1−P2A2where Piis the pressure in the chamber iand Aiis the effective section of the piston. It can be expressed as: F=P1πD2 1 4−P2π(D2 1−D2 2) 4(4.1) The force performed by the cylinder in steady-state conditions depends on whether the movement is done forward or backward, since the section is different. Assuming P2=Pr= 0 and P1=Ps, the forward force can be expressed as : Ff=PsπD2 1/4 (4.2) Concerning the backward force, P2=Psand P1=Pr= 0. The force yields: Fb=PsπD2 1−D2 2/4 (4.3) Assuming quasistatic behaviour the work can be obtained assuming the force is constant during the time and therefore multiplying the force times the displacement. 4.3 Step 3. Limiting quantities. The maximum allowed shear stress is the main quantity limiting the available mechanical force and work. It can be expressed using the Mohr circle as half the difference between the radial and tangential stresses. The radial stress in a thick walled cylinder can be written from [49] as a function of the position rin the wall as follows: σrr =PD2 1 D2−D2 11−D2 4r2(4.4) 40 4.3 Step 3. Limiting quantities. The tangential stress in a thick walled cylinder from [49] yields: σθθ =PD2 1 D2−D2 11 + D2 4r2(4.5) The equivalent shear stress can be derived from (4.4) and (4.5) as: τeq =σrr −σθθ 2=PD2 1 D2−D2 1 D2 4r2(4.6) It can be clearly seen in (4.6) that the maximum shear stress is produced for the minimum value of r, i.e. r=D1/2. Using the defined geometric relationships the maximum shear stress yields: τeq =P 1−k2 D1 (4.7) Hence, to not overcome the shear stress threshold, the maximum pressure must be established as: PL1=τeq 1−k2 D1(4.8) For backward motion, there arises another fact: there exists a maximum axial stress σaa in the rod attaching the load. It implies another pressure limitation: PL2=σaak2 D2(4.9) Then, the maximum pressure for backward motion PLb can be written as: PLb =min{PL1, PL2}=min{τeq 1−k2 D1, σaak2 D2}(4.10) Defining ϕ=σaa/τeq, ϕ > 0, it may be expressed as: PLb =τeqmin{1−k2 D1, ϕk2 D2}(4.11) 41 4. Application to hydraulic actuators 4.4 Step 4. Maximum force, stroke and work. 4.4.1 Forward motion Using (4.8) and (4.2), the maximum available force per cross-section can be expressed for the forward motion as: Ff πD2/4=τeq 1−k2 D1k2 D1(4.12) The design factor qfcan be defined as: qf=1−k2 D1k2 D1(4.13) and is the factor to be maximized in the design. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.05 0.1 0.15 0.2 0.25 kD1 = 2−1/2 Factor kD1 Design factor qf kD1 = 2−1/2 Figure 4.3: Forward force design factor Analyzing the latter expression, it can be seen that for a given size the forward force is maximized for kD1= 1/√2 performing a force per cross 42 4.4 Step 4. Maximum force, stroke and work. section of τeq/4 with a design factor qf= 1/4. Graphical results may be seen in Fig. 4.3. 4.4.2 Backward motion Using (4.11) and (4.3), the maximum available force per cross-section can be expressed for the backward motion as: Ff πD2/4=τeqmin{1−k2 D1, ϕk2 D2}k2 D1−k2 D2(4.14) An alternative formulation yields: Ff πD2/4=(τeq (1 −k2 D1) (k2 D1−k2 D2) 1 −k2 D1< ϕk2 D2 τeqϕk2 D2(k2 D1−k2 D2) 1 −k2 D1≥ϕk2 D2 (4.15) The design factor qf= 4Ff/τeqπD2may be defined as: qf=((1 −k2 D1) (k2 D1−k2 D2) 1 −k2 D1< ϕk2 D2 ϕk2 D2(k2 D1−k2 D2) 1 −k2 D1≥ϕk2 D2 (4.16) Analyzing the expression (4.16), the maximum design factor may be found by using 1 −k2 D1=ϕk2 D2or its equivalent formulation kD2=p(1 −k2 D1)/ϕ. In such a case: qf=1−k2 D1k2 D1−1 + k2 D1/ϕ(4.17) It can be expressed as: qf=−ϕ+ 1 ϕk4 D1+2ϕ+ 1 ϕk2 D1−1 (4.18) To obtain the maximum design factor: ˙qf=−4ϕ+ 1 ϕk3 D1+ 22ϕ+ 1 ϕkD1→4ϕ+ 1 ϕk3 D1max = 22ϕ+ 1 ϕkD1max (4.19) 43 4. Application to hydraulic actuators The regression analysis of each class of actuator shows that the relationship between the force and the cross-section can be found as the considered operating pressure. 102103104 100 101 102 103 Cross Section [mm2] Output force [N] Real values P=210 bar Real values P=680 bar Real values P=200 bar Maximum Optimum for steel G10100 HR Figure 4.8: Industrial hydraulic actuator force-area performance 4.8 Conclusions The present chapter has dealt with the optimization of linear hydraulic actuators. The procedure presented in [14] has been employed to obtain the maximum energy and force in a given volume, weight or cross-section. The scalability of the analyzed actuators has been also discussed. The limit output quantities have been found for certain aspect ratios and geometric relationships. The results have been compared with the performance of industrial actuators and it has been noted that the industrial actuators behave as expected. Regressions linking the force and work with 50 4.8 Conclusions the cross-section and volume have been carried out, resulting in similar performances as a function of the size as the theoretically developed models. 51 4. Application to hydraulic actuators 52 Chapter 5 Conclusions The present part has focused on providing the detailed analysis of different actuators using a general procedure and oriented towards improving the actuator design. It has introduced a new methodology to analyze linear electromagnetical and hydraulic actuators by modeling their maximum output mechanical quantities (force, work and stroke) as functions of the geometry and material properties and has discussed the scalability (in the sense of producing the same stress and strain distribution for different sizes). The motivation to undertake such a work stems from the need for light and volume reduced structures and systems, which are to be integrated in the design procedure as early as possible. Hence, the geometric relationships, aspect ratios and material properties that maximize the actuator output quantities with a certain limited volume or weight, along with their scalability for the integration in structures have been studied. A validation of the results has been done by performing dimensional analysis of the expressions obtained and comparing numerical results with industrial actuator data. 5.1 Contributions The main contribution of thesis part is the methodology described in Chapter 2along with its application to electromagnetic (Chapter 3) and hydraulic (Chapter 4) actuators. Furthermore, the contributions may be summarized as: 53 5. Conclusions •Design of a methodology to deal with the modeling and optimization of industrial actuators •Application of the methodology to linear electromagnetic actuators •Application of the methodology to linear hydraulic actuators •Design Optimization an scalability Analysis of the analyzed actuators. •Validation with real actuators and with dimensional analysis of the analyzed actuators. Part of these contributions are collected in [14] and [15]. 5.2 Future work The optimization methodology has been performed for linear hydraulic and electromagnetic actuators considering static behavior. Further research beyond the scope of this thesis is encouraged. It may be particularly interesting to investigate in the following lines: •Application of the methodology to other classes of actuators, including some other classical actuators (pneumatic) and new smart actuators (piezoelectric, magnetostrictive, shape memory alloys, magnetorheological, etc.) •Discussion of the validity of the methodology for non-linear motion. Application to rotative motion, substituting the force by the torque and the position by the angle. •Expand the methodology to deal with dynamics, taking into account other quantities, such as the power or the speed. 54 Part II Identification and Control of Piezoelectric Actuators 55 Chapter 6 Introduction Although the so-called classical actuators (electromagnetic, hydraulic and pneumatic) are the most used in the industry, new technologies based on different physical principles are being developed. In applications where the size of the actuator has to be minimized [14], or where fast response and high resolution are needed, the classical actuators fail to respond appropriately. For this reason, non-classical technologies are becoming more relevant. Among them, the piezoelectric actuators are proving to be a reliable solution for many engineering applications, ranging from micropositioning (machine tools, optic devices or modern microscopes) to active control of structures. The piezoelectric actuators are based on the known piezoelectric effect described in 1880 by Jacques and Pierre Curie [8]: in certain materials with crystalline non-symmetrical structure, dipoles are formed when the material is deformed, i.e. a mechanical strain produces an electrical field; reciprocally, the application of an electric field produces a strain. These actuators show a fast reaction time, a high resolution, a high energy density and an easy miniaturization. However, the piezoelectric actuators have some drawbacks: the reduced strain (<0.2 %), the presence of non-linearities and the high voltage needed for optimal performance. In this thesis part, we focus on the nonlinear behavior of piezoelectric actuators by taking into account the presence of hysteresis. In materials, the hysteresis is referred to the memory nature of inelastic systems where the restoring force depends not only on the instantaneous deformation but also on the history of that deformation. 57 6. Introduction In this frame, the main motivation is to delve into models to represent the hysteretic behavior of piezoelectric actuators in order to apply them to the conception of controllers for such hysteretic systems. These controllers can allow a more optimum control of the devices employing piezoelectric actuators. 6.1 Modeling and validation of piezoelectric actuators It is known that the presence of non-linearities and the high voltage needed for optimal performance are the main drawbacks of piezoelectric actuators. We focus on the nonlinear behavior of piezoelectric actuators by taking into account the presence of hysteresis. To describe the behavior of hysteretic processes several mathematical models have been proposed [54]: the Duhem model [9] uses the property that a hysteretic system’s otput changes its character when the input changes direction; the Ishlinskii hysteresis operator has been proposed as a model for plasticity-elasticity [36]; the Preisach model has been used for the modeling of electromagnetic hysteresis [43]; the Bouc-Wen model has been used to model wood joints and structural systems [11]. A survey of the mathematical models for hysteresis may be found in [40]. These models have been applied to describe the behavior of piezoelectric actuators: Prandtl-Ishlinskii in [44], Preisach in [51] and Bouc-Wen in [39]. An energy based model has been employed in [45]. The present is focused on the Bouc-Wen model. In a recent work, the hysteresis loop obtained by the Bouc-Wen model has been characterized analytically [29]. In further work [26], a new parametric nonlinear identification technique for the Bouc-Wen model based on the analytical description of [29] is proposed. This method does not use any information from the behavior of the system in the plastic region which makes it applicable for a wide class of materials including base isolation devices, magnetorheological 58 6.2 Control of piezoelectric actuators considering the hysteresis damper, piezoelectric elements, etc. And, unlike most identification techniques for the Bouc-Wen model, this method provides the exact values of the model parameters in the absence of disturbances, and gives a guaranteed relative error between the estimated parameter and the true ones in the presence of a particular class of perturbations. The main advantages of the proposed identification methodology are (1) the simplicity of the proof that the estimated parameters are within a given tolerance with respect to their true counterparts in the presence of disturbances; (2) the fact that limit cycles can be obtained experimentally in a simple way [4]; (3) its wider range of applicability than [30]. The identification technique consists in exciting the hysteretic system with two periodic signals that have a specific shape. The parameters of the Bouc-Wen model are then obtained from the two limit cycles using a precise algorithm. The method guarantees that the estimated parameters are within a given tolerance with respect to the true parameters, and it is shown that the identification technique is robust with respect a class of disturbances of practical interest. However, when applying the method to the modeling of certain piezoelectric actuators, an inexact matching has been noted. To improve the matching between the model and the experimental behaviour of a certain piezoelectric actuator, we propose a modification of the Bouc-Wen model. To identify such a modified model, we have developed a new identification technique based on the results obtained in [26]. The modified Bouc-Wen model is validated by means of experiments, and compared to the behavior of the non-modified Bouc-Wen model. 6.2 Control of piezoelectric actuators considering the hysteresis The main challenge for the control design of applications with piezoelectric actuators is the presence of hysteresis. In this work, we consider the problem of micropositionning using a piezoelectric actuator. This problem has 59 7. Piezoelectricity It can be expressed in a single expression as:               S1 S2 S3 S4 S5 S6 D1 D2 D3               =               sE 11 sE 12 sE 13 0 0 0 0 0 d31 sE 12 sE 11 sE 13 0 0 0 0 0 d31 sE 13 sE 13 sE 33 0 0 0 0 0 d33 000sE 44 0 0 0 d15 0 0000sE 44 0d15 0 0 00000sE 66 000 0000d15 0ε11 0 0 000d15 0 0 0 ε11 0 d31 d31 d33 0 0 0 0 0 ε33               ×               T1 T2 T3 T4 T5 T6 E1 E2 E3               (7.9) The so-called electromechanical coupling factor kis specially significative in the characterization of a piezoelectric element, it is defined in [52] as: k2=d2 sEεT(7.10) and it shows the relationship between the stored mechanical energy and the input electrical energy when working as an actuator, and between the stored electrical energy and the input mechanical energy when working as a sensor. 7.1.1 A brief history Some references [19;35] deal with the history of piezoelectric technology. The most important events are reported here briefly. In 1880, the brothers Pierre Curie and Jacques Curie predicted and demonstrated piezoelectricity using tinfoil, glue, wire, magnets, and a jeweler saw [8;58]. They showed that crystals of tourmaline, quartz, topaz, cane sugar, and Rochelle salt (sodium potassium tartrate tetrahydrate) generate electrical polarization from mechanical stress. Quartz and Rochelle salt exhibited the most piezoelectricity. The term piezoelectricity was first suggested by W. Hankel in 1881. Converse piezoelectricity was mathematically deduced from fundamental thermodynamic principles by Lippmann in 1881. The Curies immediately confirmed the existence of the converse effect and obtained quantitative proof of the complete reversibility of deformations in piezoelectric crystals. 66 7.1 The piezoelectric effect In 1910 Voigt published Lehrbuch der Kristallphysik [55], and it became a standard reference work detailing the complex electromechanical relationships in piezoelectric crystals. During World War I the piezoelectric ultrasonic transducer was developed by Langevin. At the same time piezoelectric materials began to be used as microphones, accelerometers, underwater transducers, etc. However, the limited material performance inhibited commercialization. During World War II BaTiO3was discovered as a high dielectric constant material in USA, UK, USSR, and Japan, independently. Gray discovered a poling process, which made ceramic materials act as a single crystal possessing both ferroelectric and piezoelectric properties. In 1952, PZT was reported as ferroelectric solid-solution system, and the phase diagram was established by Shirane, et al. PZT was reported as useful piezoelectric transducer material by B. Jaffe et al. in 1954. Piezoelectric ceramics applications became commercialized, including phonograph pick-ups, microphones, underwater transducers (sonar), ignition systems, discrete actuators, etc. The 1960s-1980s decades were important for the discovery and research of transparent electro-optic (Pb, La)(Zr, Ti)O3PLZT ceramics and by the development of Pb(Mg1/3Nb2/3)O3PMN and other relaxor ferroelectric ceramics and devices. Also there was the first development of multi-layer stack actuators. From 1980 to now piezoelectric actuators has been used for smart structures, distributed actuator systems, prototype smart beam, active airfoil, etc.. There has been a development of flexible actuators based on piezoelectric fibers embedded in polymer matrix (active fiber composites), with applications for active vibration reduction and noise control system. The use of piezoelectric actuators in micro and nano positioning devices requiring high precision such as modern microscopes is one of the actual applications and challenges of the piezoelectric technology. 67 7. Piezoelectricity 7.1.2 Deformation modes. The deformation directions are shown in Fig. 7.1. It is important to note that all the parameters used in (7.1) have to be considered in the different deformation directions. Figure 7.1: Axes and deformation directions. Depending on the electrical field application and the deformation of interest, piezoelectric actuators can be employed using different modes: •Longitudinal mode d33. See Fig. 7.2(a). Expression (7.6) turns into: S3= 6 X i=1 sE 3iTi+d33E3(7.11) •Transverse mode d31. See Fig. 7.2(b). Expression (7.6) turns into: S1= 6 X i=1 sE 1iTi+d31E3(7.12) •Shear mode d15. See Fig. 7.2(c). Expression (7.6) turns into: S5= 6 X i=1 sE 5iTi+d15E1(7.13) 68 7.2 Piezoelectric actuator simplified model (a) Longitudinal mode. (b) Transverse mode. (c) Shear mode. Figure 7.2: Different deformation modes. 7.2 Piezoelectric actuator simplified model 7.2.1 Low frequency A piezoelectric element can be modeled from (7.1) as the association in parallel of a capacitor and a charge source, since the charge can be obtained from the electric displacement D, and the voltage can be derived from the electrical field E, assuming that it is uniformly distributed in a length l(V=E/l). Expression (7.1) can be written as: Qe A=dF A+εTV z x l0 =sEF A+dV z (7.14) where Qeis the electrical charge, Ais the cross-section in the movement direction, Fis the force, Vis the applied voltage, l0is the initial length in the movement axis and zis the thickness in electrical field direction. The first equation of (7.14), (known as the sensor expression) can be written as: Qe=dF +εTA zV=dF +CV (7.15) where C=εTA/z is the equivalent capacitance. The second equation of (7.14), (known as the actuation expression) can be written as: x=sEl0 F A+V dl0 z=k−1F+dl0 zV(7.16) where k=A/sEl0is the equivalent stiffness constant. Note that in the longitudinal mode, the electrical field is applied in the motion’s direction 69 7. Piezoelectricity and thus: x=k−1F+dV (7.17) 7.2.2 High frequency These approximations of (7.15),(7.16) and (7.17) apply for low frequencies but when the dynamic behavior for higher frequencies (close to the mechanical resonance frequency) is concerned, the model from [56] characterized in Fig. 7.3 has to be used. It includes the equivalent capacitor and a RLC branch in parallel where R1includes the mechanical losses, L1is the equivalent inductance of the mechanical circuit and C1the capacitance of the mechanical circuit. Each branch has a mechanical resonance at fi= 1/2π√LiCi. A current (or charge) source can be added if the system is mechanically loaded. More branches can be added corresponding to the resonance frequencies of the mechanical system. Figure 7.3: Equivalent circuit of a piezoelectric element excited at high frequency. The impedance behavior against the frequency considering only one resonance frequency is plotted in Fig. 7.4. It can be noted that the frequencies between them show an inductive behavior while the others below resonance and above antiresonance show capacitive behavior. The resonance and antiresonance frequencies can be found for values close to the series and parallel 70 7.2 Piezoelectric actuator simplified model 102103104105 10−1 100 101 102 103 104 Frequency (Hz) Impedamce (Ω) R=0.1Ω R=1Ω R=10Ω (a) Impendance value 102103104105 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 Frequency (Hz) Impedamce angle(rad) R=0.1Ω R=1Ω R=10Ω (b) Impedance angle Figure 7.4: Impedance of a piezoelectric element with different R values and C0= 0.1µF, C1= 1 µF and L1= 0.1 mH. 71 7. Piezoelectricity resonant frequency as: fr=1 2πr1 L1C1 fa=1 2πrC0+C1 L1C0C1 (7.18) 7.2.3 Load The relationship between force and displacement can be extracted from expression (7.16). Manufacturers usually provide the force with no displacement and the free displacement. Defining F0as the force with no displacement (clamped actuator) and x0the free displacement with no force: x0=dV l0 z F0=dV A zsE (7.19) Hence expression (7.16) can be rewritten as: F=F0 x0 (x0−x) (7.20) where both F0and x0depend linearly on the applied voltage. Note that the previously defined stiffness constant k, can be expressed as F0/x0and does not depend on the voltage but on the material stiffness. An alternative expression of (7.20) is: F=k(x0−x) = F0−kx (7.21) 7.2.3.1 Example An example can be shown with a sample actuator working in the transversal mode. The parameters are: 72 7.2 Piezoelectric actuator simplified model l0= 50 ·10−3m z0= 0.2·10−3m A= 6 ·10−6m2 sE 31 = 15 ·10−12m2/N d31 =−250 ·10−12m/V Then: k=A sE·l0 =6·10−6 15 ·10−12 ·50 ·10−3= 8 ·106N/m For V= 400 V: x0=d·V·l0 z0=−250 ·10−12 ·400 ·50·10−3 0.2·10−3= 25 ·10−6m F0=x0·k= 200N 0 0.5 1 1.5 2 2.5 x 10−5 0 20 40 60 80 100 120 140 160 180 200 Displacement [m] Force [N] V= 80V V= 160V V= 240V V= 320V V= 400V Constant load Elastic load Figure 7.5: Displacement - Force curves In Figure 7.5 the load - displacement characteristic for different voltages can be seen. Also the load - displacement characteristic for different voltages under a constant load and linear load (for example a spring or a attached structure) are shown. 73 7. Piezoelectricity 7.3 Considerations 7.3.1 Non-linearities It is known that the presence of non-linearities is one of the main drawbacks of piezoelectric actuators. The most important non-linearities involved in piezoelectric materials are hysteresis and creep. The hysteresis is referred to the memory nature of inelastic systems where the restoring force depends not only on the instantaneous deformation but also on the history of that deformation. The hysteresis (Fig. 7.6) is produced by the retarded reorientation of dipole domains, which initially maintain their direction in the field direction upon reducing its strength. The creep refers to the time variation of the strain. When the field is changed and hold constant at a certain level, more and more dipoles orient themselves in the applied direction and thus, a increment in the strain is produced. It is important to note that creep is significative in static conditions while hysteresis has to be always taken into account. To describe the behavior of hysteretic processes several mathematical models have been proposed [54]: the Duhem model [9] uses the property that a hysteretic system’s otput changes its character when the input changes direction; the Ishlinskii hysteresis operator has been proposed as a model for plasticity-elasticity [36]; the Preisach model has been used for the modeling of electromagnetic hysteresis [43]; the Bouc-Wen model has been used to model wood joints and structural systems [11]. A survey of the mathematical models for hysteresis may be found in [40]. These models have been applied to describe the behavior of piezoelectric actuators: Prandtl-Ishlinskii in [44], Preisach in [51] and Bouc-Wen in [39]. An energy based model has been employed in [45]. In the present thesis part we consider the modeling of a piezoelectric actuator using the Bouc-Wen model for smooth hysteresis [57]. This model has received an increasing interest due to its ability to capture in an analytical form a range of shapes of hysteretic cycles which match the behavior of a wide class of hysteretical systems [46]. In particular, it has been used to 74 7.3 Considerations Figure 7.6: Example of a displacement - voltage hysteresis curve model piezoelectric elements [39], magnetorheological dampers [7;47] and wood joints [11]. The models, derived from experiments, have been used either to predict the behavior of the physical hysteretic element [47] or for control purposes as in [6;28;31]. 7.3.2 Temperature dependance The temperature is an important quantity to be considered when dealing with piezoelectric actuators. The Curie temperature TCis a threshold value. Above TCthe piezoelectric materials lose their piezoelectric properties. The Curie temperature ranges from 160 ◦Cto 350 ◦Cdepending on the materials. It is important to remark that depolarization begins to occur below TCand thus the temperature should be limited to half of the Curie temperature. 75 8. The Bouc-Wen model which is the solution of the nonlinear first order differential equation (8.2). In this equation, A, β and γare nondimensional parameters which control the shape and the size of the hysteresis loop, while nis a scalar that governs the smoothness of the transition from elastic to plastic response. The Bouc- Wen model has received an increasing interest due to its ability to capture in an analytical form a range of shapes of hysteretic cycles which match the behavior of a wide class of hysteretical systems [46]. In particular, it has been used to model piezoelectric elements [39], magnetorheological dampers [7;47] and wood joints [11]. The models, derived from experiments, have been used either to predict the behavior of the physical hysteretic element [47] or for control purposes as in [6;28;31]. 8.2 The normalized Bouc-Wen model 8.2.1 Classification of the Bouc-Wen models The nonlinear hysteretic behavior may be conceptualized [32] as a map x(t)7→ Φs(x)(t), where x(t) represents the time history of an input variable and Φs(x)(t) describes the time history of the hysteretic output variable. Two fundamental properties are shared by many physical hysteretic systems arising from structural, mechanical and electromechanical engineering: Property 1: For any bounded input x(t), the output of the true hysteresis Φs(x)(t) is bounded. This bounded input-bounded output (BIBO) property stems from the fact that, in practice, many (electro)mechanical and structural systems are stable in open loop. Property 2: The physical systems that include hysteretic components dissipate energy. To represent adequately the true hysteresis Φs(x)(t), the Bouc-Wen model ΦBW (x)(t) needs to keep both properties, that is to be BIBO and dissipative. Define the sets: 82 8.2 The normalized Bouc-Wen model Ωα,k,D,A,β,γ,n ={z(0) ∈Rsuch that ΦBW is BIBO with fixed values of the parameters α, k, D, A, β, γ, n}(8.3) ΩA,β,γ,n ={z(0) ∈Rsuch that z(t) is bounded for any C1bounded input signal x(t)with fixed values of the parameters A, β, γ, n}(8.4) Ω? A,β,γ,n ={z(0) ∈Rsuch that z(t) is bounded for any C1input signal x(t) with fixed values of the parameters A, β, γ, n}(8.5) Then, we have the following result which characterizes the two classes of Bouc-Wen models that are BIBO and asymptotically dissipative [29]. Theorem 1 Define the constants: z0,n sA β+γand z1,n sA γ−β.(8.6) Then, Table 8.1 holds. Table 8.1: Classification of the BIBO, passive and thermodynamically consistent Bouc-Wen models CASE ΩA,β,γ,n Upper bound on |z(t)|CLASS A > 0β+γ > 0 and β−γ≥0Rmax (|z(0)|, z0) I Furthermore, we have Ωα,k,D,A,β,γ,n = Ω? A,β,γ,n = ΩA,β,γ,n (8.7) A by-product of Theorem 1 is the existence and uniqueness of the solution z(t) over t∈[0,+∞). Equality (8.7) means that the boundedness of 83 8. The Bouc-Wen model the signal z(t) depends only on the parameters A,γ,βand n, and it is independent of the boundedness of the input signal x(t). This fact is particularly important for system control theory: when x(t) is a closed loop signal, we cannot assume a priori that it is bounded. The fact that Ω? A,β,γ,n = ΩA,β,γ,n shows that for every input signal x(t) (under the only assumption that it is C1), the output z(t) is always bounded if the set Ω is non-empty, and if z(0) ∈Ω. In parallel work [10], the study of the thermodynamic admissibility of the Bouc-Wen model within the context of the endochronic theory led to the following result: the conditions A > 0 and −β⩽γ⩽βare necessary and sufficient for the thermodynamic admissibility of the Bouc-Wen model. This means that the class I Bouc-Wen model is consistent with the laws of thermodynamics. 8.2.2 The normalized Bouc-Wen model Consider two Bouc-Wen models (8.1)-(8.2) whose parameters are such that n2=n1=n,A2=A1,β2=νnβ1,γ2=νnγ1,D2=νD1,α2=α1,k2=k1 where νis a positive constant, and with an initial condition z2(0) = z1(0) = 0. Then both models belong to the same class, and for any input signal x(t) they deliver exactly the same output ΦBW (t). This means that the inputoutput behavior of a Bouc-Wen model is not described by a unique set of parameters {α, k, D, A, β, γ, n}and, for this reason, identification procedures that use input-output data cannot determine the parameters of the Bouc-Wen model. To cope with this problem, users of the Bouc-Wen model often fix some parameters to arbitrary values as in reference [41] where the coefficient (1 −α)Dk of z(t) in equation (8.1) has been set to one and the parameter Dhas also been set to one. Other authors compare the shape of the limit cycle instead of comparing the identified parameters with their true values as in reference [47]. This fact makes it very difficult to compare results of different identification methods by comparing the identified parameters. Thus it is necessary to elaborate some equivalent “normalized” model whose parameters define in a unique way the input-output behavior of the model 84 8.2 The normalized Bouc-Wen model allowing a parametric-based comparison of identification methods for this hysteretic model. To this end, define w(t) = z(t) z0 so that the model (8.1)- (8.2) can be written as: ΦBW (x)(t) = κxx(t) + κww(t),(8.8) ˙w(t) = ρ˙x−σ|˙x(t)||w(t)|n−1w(t)+(σ−1) ˙x(t)|w(t)|n(8.9) where ρ=A Dz0 >0, σ =β β+γ≥0, κx=αk > 0, κw= (1 −α)Dkz0>0. (8.10) We call equations (8.8)-(8.9) the normalized form of the Bouc-Wen model. Note that if the initial condition w(0) is such that |w(0)| ≤ 1 then, by Theorem 1,|w(t)| ≤ 1 for all t≥0. This means that the variable z(t) has been scaled to unity. It can be checked that the normalized form of the Bouc- Wen model defines a bijective relationship between the input-output behavior of the model and its parameters. It also has the advantage of having only five parameters to identify instead of the seven parameters for the standard form. Note that the normalized form of the Bouc-Wen model is exactly equivalent to its standard form. Indeed, for any input x(t), both forms deliver exactly the same output ΦBW (t) taking into account that we have w(0) = z(0) z0 . The classification of the normalized Bouc-Wen models is given in Table 8.2. It Table 8.2: Classification of the BIBO, passive and thermodynamically stable normalized Bouc-Wen models CASE Ωσ,n Upper bound on |w(t)|CLASS σ≥1 2Rmax (|w(0)|,1) I can be seen that a single parameter σis needed for this classification. 85 8. The Bouc-Wen model With these notations we obtain from equation (8.9): For w(t)≥0,˙x(t)≥0 ˙w(t) = ρ(1 −w(t)n) ˙x(t) (8.11) For w(t)≤0,˙x(t)≥0 ˙w(t) = ρ(1 + (2σ−1) (−w(t))n) ˙x(t)(8.12) For w(t)≥0,˙x(t)≤0 ˙w(t) = ρ(1 + (2σ−1)w(t)n) ˙x(t) (8.13) For w(t)≤0,˙x(t)≤0 ˙w(t) = ρ(1 −(−w(t))n) ˙x(t) (8.14) 86 Chapter 9 Analysis and parameter identification of the Bouc-Wen model It is known that the presence of non-linearities and the high voltage needed for optimal performance are the main drawbacks of piezoelectric actuators. As it is explained in Section 7.3.1, we focus on the nonlinear behavior of piezoelectric actuators by taking into account the presence of hysteresis. The normalized version of the model introduced in chapter 8relates the output restoring force ΦBW (x)(t) to the input displacement x(t) in the following way: ΦBW (x)(t) = κxx(t) + κww(t),(9.1) ˙w(t) = ρ˙x(t)−σ|˙x(t)||w(t)|n−1w(t)+ +(σ−1) ˙x(t)|w(t)|n) (9.2) where κx>0, κw>0, ρ > 0, σ > 1 2and n≥1 are the model parameters that shape the hysteresis loop. The range of the parameter σcorresponds to the Class I Bouc-Wen model which is stable, asymptotically dissipative and thermodynamically consistent [29]. This chapter deals with the problem of identifying the model parameters in the presence of disturbances. The signals that are accessible to measurements are the input x(t) and the output 87 9. Analysis and parameter identification of the Bouc-Wen model ΦBW (x)(t). The state w(t) is not accessible to measurements. As can be seen from equations (9.1)-(9.2), the difficulty of the identification problem lies (1) in the nonlinear form of the model, especially in relation with the estimation of the parameter nwhich forms part of the “structure” of the model and (2) in the fact that the state w(t) is not accessible to measurements. A survey of the parametric and non parametric methods that have been used in the literature for the identification of the Bouc-Wen model may be found in [38]. The main theoretical deficiency of these methods is that they rely mainly on numerical simulations and do not offer, to a large extent, a rigorous mathematical proof of the convergence of the estimated parameters to their true counterparts. In this chapter, we propose a new parametric nonlinear identification technique for the Bouc-Wen model based on the analytical description of [29]. This method does not use any information from the behavior of the system in the plastic region which makes it applicable for a wide class of materials including base isolation devices, magnetorheological damper, piezoelectric elements, etc. And, unlike most identification techniques for the Bouc-Wen model, this method provides the exact values of the model parameters in the absence of disturbances, and gives a guaranteed relative error between the estimated parameter and the true ones in the presence of a particular class of perturbations. The main advantages of the proposed identification methodology are (1) the simplicity of the proof that the estimated parameters are within a given tolerance with respect to their true counterparts in the presence of disturbances (2) the fact that limit cycles can be obtained experimentally in a simple way [4] (3) its wider range of applicability than [30]. The identification technique consists in exciting the hysteretic system with two periodic signals that have a specific shape. The parameters of the Bouc-Wen model are then obtained from the two limit cycles using a precise algorithm. This method guarantees that the estimated parameters are within a given tolerance with respect to the true parameters, and it is shown that the identification technique is robust with respect a class of disturbances of practical interest. 88 9.1 Parameter identification for the Bouc-Wen model 9.1 Parameter identification for the Bouc-Wen model 9.1.1 Class of inputs In this chapter we consider that the input signal x(t) is T-wave periodic [29]. This means that it is continuous on the time interval [0,+∞) and periodic of period T > 0. Furthermore there exists a scalar 0 < T+< T such that the signal xis C1on both intervals (0, T+) and (T+, T) with ˙x(τ) = dx(τ) dτ >0 for τ∈(0, T+) and ˙x(τ)<0 for τ∈(T+, T) (see Figure 9.1). We denote Xmin =x(0) and Xmax =x(T+)> Xmin the minimal and maximal values of the input signal, respectively. We assume that max (|Xmax|,|Xmin|)≤κw κx so that the Bouc-Wen model is consistent with the hysteretic property [30]. Figure 9.1: Example of a T-wave periodic signal. 89 9. Analysis and parameter identification of the Bouc-Wen model 9.1.2 Analytic description of the forced limit cycle for the Bouc-Wen model Define the following functions: ϕ− σ,n(w) = Zw 0 1 1 + σ|u|n−1u+ (σ−1)|u|ndu (9.3) ϕ+ σ,n(w) = Zw 0 1 1−σ|u|n−1u+ (σ−1)|u|ndu (9.4) ϕσ,n(w) = ϕ+ σ,n(w) + ϕ− σ,n(w) (9.5) for any scalar w∈(−1,1). In this section and in the rest of the chapter we denote w(t) the solution of the differential equation (9.2) while the notation wwithout an argument is used for a given scalar. It has been shown in [29] that the functions ϕ− σ,n(·), ϕ+ σ,n(·) and ϕσ,n(·) are strictly increasing on the interval (−1,1) so that they are bijective. Their inverses are denoted ψ− σ,n(·), ψ+ σ,n(·) and ψσ,n(·), respectively. These functions have been studied extensively in [29]. Note that for w≥0 we have ϕ− σ,n(w) = Zw 0 1 1 + (2σ−1)undu (9.6) ϕ+ σ,n(w) = Zw 0 1 1−undu (9.7) and for w≤0 we have ϕ− σ,n(w) = Zw 0 1 1−(−u)ndu (9.8) ϕ+ σ,n(w) = Zw 0 1 1 + (2σ−1)(−u)ndu (9.9) The limit cycle for the Bouc-Wen model is described by the following [29]: Theorem 2 Let x(t)be a T-wave periodic input signal. Define the functions ωmand φmfor any positive integer mas follows ωm(τ) = w(mT +τ)for τ∈[0, T] (9.10) φm(τ) = κxx(τ) + κwωm(τ)for τ∈[0, T] (9.11) 90 9.1 Parameter identification for the Bouc-Wen model where w(·)is the solution of equation (9.2) with initial condition w(0). Then the sequence of functions {φm}m≥1(resp. {ωm}m≥1) converges uniformly on the interval [0, T]to a continuous function ¯ ΦBW (resp. ¯w) defined as ¯ ΦBW (τ) = κxx(τ) + κw¯w(τ)for τ∈[0, T] (9.12) ¯w(τ) = ψ+ σ,n ϕ+ σ,n [−ψσ,n (ρ(Xmax −Xmin))] + +ρ(x(τ)−Xmin)) for τ∈[0, T+] (9.13) ¯w(τ) = −ψ+ σ,n ϕ+ σ,n [−ψσ,n (ρ(Xmax −Xmin))] −ρ(x(τ)−Xmax)) for τ∈[T+, T] (9.14) Furthermore we have for all τ∈[0, T] −1<−ψσ,n (ρ(Xmax −Xmin)) ≤¯w(τ) ≤ψσ,n (ρ(Xmax −Xmin)) <1 (9.15) the lower and upper bounds of ¯w(τ)being attained at τ= 0 and τ=T+ respectively. 9.1.3 Identification methodology In general, the nonlinear state variable wis not accessible to measurement. However, in many cases of practical importance, the hysteretic limit cycle can be obtained experimentally [41]. The hysteretic system under study is assumed to be described by the normalized Bouc-Wen model (9.1)-(9.2), with unknown parameters κx,κw,ρ,σand n. The loading part of the limit cycle (that corresponds to an increasing input x(t)) can be obtained from Theorem 2as: ¯ ΦBW (x) = κxx+κw¯w(x) (9.16) ¯w(x) = ψ+ σ,n ϕ+ σ,n [−ψσ,n (ρ(Xmax −Xmin))] + +ρ(x−Xmin)) (9.17) 91 9. Analysis and parameter identification of the Bouc-Wen model tions, it is also desirable that a “small size” of disturbances leads to a “small discrepancy” between the identified parameters and the true one. Theorem 3 says that, given an  > 0, for all µ-small disturbances such that 0 ≤µ≤µ∗, the relative error between the identified parameters p◦and the true parameters pdoes not exceed . If the quantity µ∗were zero, this would have implied that, even for arbitrarily small disturbances, the identification method may lead to a large discrepancy between the identified parameters and the true ones. Theorem 3guarantees that the robustness margin µ∗>0 so that all µsmall disturbances with µ∈[0, µ∗] lead to a relative error in the parameters no more than . 9.2 Numerical simulation example In this section we consider the Bouc-Wen model given by the unknown parameters κx= 2, κw= 2, ρ= 1, σ= 3, n= 1.5. The objective is to use the technique presented in the previous sections to identify its parameters. As seen in Section 9.1.3, the identification technique has 10 steps. Step 1. The first step of the identification procedure is the choice of the T-periodic input signals. Due to Assumption 1we have  ˙ ξ(τ)≤µ|˙x(τ)|. This implies that derivative ˙ ξ(τ) of the disturbance ξ(τ) needs to be zero whenever the derivative of the input signal x(τ) is zero. Thus, a sine wave input signal candidate would impose that ˙ ξ(τ) should be very small around the time instants 0 + mT and T 2+mT (mis any positive integer) which is unlikely to happen in practice. For this reason, a good choice of an input signal is a triangular one so that the derivative ˙ ξ(τ) needs only to be small with respect to the slope of the input signal which is constant (in absolute value). The next design parameter to be chosen is the frequency of the input signal. Since the Bouc-Wen model is rate independent, its input-output behavior is independent of the frequency of the input signal. We thus take T= 1 and T+=T 2. We also choose Xmax =−Xmin = 0.2. Step 2. In this step, one has to choose a value q6= 0 to obtain a second input signal x1(t) = x(t) + q. The signals x(t) and x1(t) are given in Figure 98 9.2 Numerical simulation example 9.2 with q= 0.1. 024 −0.5 0 0.5 Input signals Time 024 −1 0 1 2 Output signals Time −0.5 0 0.5 −1 −0.5 0 0.5 1 1.5 Displacement Force Figure 9.2: Upper left. Solid: input signal x(t), dashed: input signal x1(t). Lower left. Solid: output ΦBW (x)(t), dashed: output ΦBW,1(x)(t). Right. Limit cycles (x, ¯ ΦBW ) (solid) and (x1,¯ ΦBW,1) (dashed) that have been obtained for the time interval [4T, 5T] In practice, the input and output data are in the form of a finite number of samples x(kh), ¯ Φ(kh) where his the sampling period, k= 0,1,··· , m and mthe number of samples. These samples have to be taken once the output of the system is in steady-state. Note that, since the identification technique uses only the loading part of the limit cycle, we can choose the time instant kh = 0 such that x(k= 0) corresponds to the lowest value of xand the time instant mh so that x(k=m) corresponds to the largest value of x. This implies that the samples that are used for identification purpose verify 99 9. Analysis and parameter identification of the Bouc-Wen model x(i)< x(i+ 1) for all 0 ≤i < m as we are considering the loading part of the limit cycle. Step 3. The estimate κ◦ xof the coefficient κwis computed from equation (9.20) as: κ◦ x=¯ ΦBW,1(x(0) + q)−¯ ΦBW (x(0)) q(9.37) where x(0) is the value of xat the time instant k= 0. Step 4. An estimate θ◦(x) of the function θ(x) is computed from equation (9.21) as θ◦(x(i)) = ¯ ΦBW (x(i)) −κ◦ xx(i) for i= 0,··· , m (9.38) Step 5. It has been shown in the previous section that the estimate θ◦(x) is strictly increasing and has a unique zero, that is there exists a unique point x∗such that θ◦(x∗) = 0. Since all the samples x(i) are such that x(i)< x(i+ 1), we have θ◦(x(i)) < θ◦(x(i+ 1)). The existence and unicity of the zero of the function θ◦shows that there exists a unique integer rsuch that θ◦(x(r)) ≤0< θ◦(x(r+ 1)). This implies that x(r)≤x∗< x(r+ 1), and a linear interpolation gives an estimate x◦ ∗of the zero x∗. A simple computer program can be done to determine the integer r. Step 6. An estimate of the parameter ais computed from equation (9.23) as: a◦=θ◦(x(r+ 1)) −θ◦(x(r)) x(r+ 1) −x(r)(9.39) Step 7 Choosing the design parameters x∗2=x(l2)> x∗1=x(l1)> x◦ ∗, the estimates n◦and b◦are computed from equations (9.24) and (9.25) as follows: n◦= log θ◦(x(l2+1))−θ◦(x(l2)) x(l2+1)−x(l2)−a◦ θ◦(x(l1+1))−θ◦(x(l1)) x(l1+1)−x(l1)−a◦! log θ◦(x∗2) θ◦(x∗1)(9.40) b◦=a◦−θ◦(x(l2+1))−θ◦(x(l2)) x(l2+1)−x(l2) θ◦(x∗2)n◦(9.41) Step 8 Estimates of the parameters κwand ρare computed from equations 100 9.3 Conclusion (9.26) and (9.27) as follows: κ◦ w=n◦ ra◦ b◦(9.42) ρ◦=a◦ κ◦ w (9.43) Step 9 An estimate of the function ¯w(x) is computed from equation (9.28) as follows: ¯w◦(x(i)) = θ◦(x(i)) κ◦ w for i= 0,··· , m (9.44) Step 10 Choose a design parameter x∗3=x(l3)< x◦ ∗. Then an estimate of the parameter σis computed from equation (9.29) as: σ◦=1 2         ¯w◦(x(l3+ 1)) −¯w◦(x(l3)) x(l3+ 1) −x(l3) ρ◦−1 (−¯w◦(x∗3))n◦+ 1         (9.45) The numerical simulation gives κ◦ x= 2.0000, κ◦ w= 2.0059, ρ◦= 0.9971, n◦= 1.4954, σ◦= 2.9728. 9.3 Conclusion This chapter has presented a new identification method for the Bouc-Wen model. The method consist in exciting the hysteretic systems with two input signals that differ by a constant, and use the obtained limit cycles to derive the parameters of the Bouc-Wen model. This technique provides the exact values of the parameters in the absence of disturbances, and proves to be robust with respect to a class of perturbations of practical relevance. 101 9. Analysis and parameter identification of the Bouc-Wen model 102 Chapter 10 Adaptation of the Bouc-Wen model for the modeling and validation of a piezoelectric actuator In this chapter, we propose a modification of the Bouc-Wen model to describe the experimentally observed behavior of a piezoelectric actuator. To identify this modified model, we have developed a new identification technique based on the results obtained in [26], where the problem of identifying the Bouc-Wen model parameters is addressed. The modified Bouc-Wen model is validated by means of experiments, and is compared to the behavior of the non-modified Bouc-Wen model. The chapter is structured as follows. Section 10.1 shows that the model presented in the last chapter does not describe with precision the experimental behavior of the piezoelectric actuator. In Section 10.2 the modified Bouc-Wen model is introduced, along with the corresponding parameter identification methodology. Section 10.3 applies the identification technique of Section 10.2 and validates the obtained model using experiments. It also presents a comparison between the modified and non-modified Bouc-Wen models. The conclusions are summarized in section 10.4. 103 10. Adaptation of the Bouc-Wen model for the modeling and validation of a piezoelectric actuator 10.1 Experimental observations The system under study is the patch of Figure 10.1 which is a piezoelectric actuator that contains the foil PIC-255 (Physik Instrumente). The actuator is seen as a SISO system whose input is the voltage vapplied to the 3 axis and the output is the displacement yalong the 1 axis. The model of the piezoelectric actuator is given by: Figure 10.1: Piezoelectric patch employed for the experiments. m¨y(t) + c˙y(t) + k1(y(t)−y0) + k2w(t) = k3v(t) (10.1) where mis the equivalent mass of the free edge point of piezoelectric actuator, y(t) its relative position with respect to the sensor, y0is a constant that depends on the choice of the origin, vthe input voltage, and ki,i= 1,2,3 are constant gains. The nonlinear term w(t) takes into account the effect of hysteresis. We use in this section periodic input voltage functions that have a low frequency. In this case, the terms m¨y(t) and c˙y(t) may be neglected so 104 10.2 The modified model and the corresponding identification methodology that the model of the piezoelectric actuator can be written as: y(t) = kvv(t) + kww(t) + y0(10.2) where kvand kware constant gains. In the rest of the section, we approximate the nonlinear term w(t) with a Bouc-Wen model and we use the identification method of Section 9.1.3 to determine its parameters. Note that the input variable is the voltage v(which plays the role of xin equations (9.1)-(9.2)), and the output variable is y(which plays the role of Φ(x) in equations (9.1)- (9.2)). According to this methodology, two wave T-periodic voltages of low frequency f= 0.1Hz, and that differ by a constant offset of q= 100 Vare applied to the actuator. Figure 10.2 upper gives the two limit cycles obtained asymptotically as a response of the actuator to the two input voltages. If the actuator were described precisely by the Bouc-Wen model, we would have from equation (9.20): ¯ Φ1(v+q)−¯ Φ(v) = κvq(10.3) for any value of v∈[Vmin, Vmax]. Hence, such a difference would be constant so that a drag-and-drop of the two voltage-displacement curves of Figure 10.2 upper would lead to a perfect matching. However, we observe in Figure 10.2 lower that this is not the case. This means that the model composed of equations (9.1)-(9.2), (10.2) does not describe satisfactorily the experimental behavior of the piezoelectric actuator. The next section is dedicated to modifying this model so that it matches with experimental observations. 10.2 The modified model and the corresponding identification methodology In the previous section, it has been observed that the Bouc-Wen model does not represent precisely the experimental behavior of the piezoelectric actuator. For this reason, we propose a modification of the model which consists in 105 10. Adaptation of the Bouc-Wen model for the modeling and validation of a piezoelectric actuator 0 50 100 150 200 250 300 350 400 0 10 20 30 40 Voltage [V] Displacement [μ m] 0 50 100 150 200 250 300 350 400 0 10 20 30 40 Voltage [V] Displacement [μ m] Figure 10.2: Drag and drop of the voltage-displacement curve of 100−400 V input signal over the 0 −300 Vinput signal. It can be seen that the curves do not match. introducing a higher degree polynomial in the input variable instead of a linear term. We also propose a modification of the identification methodology of Section 9.1.3. 10.2.1 Modified model The term κxx(t) of (9.1) is substituted by a polynomial function as: Φ(x)(t) = N X i=1 κixi(t) + κww(t) (10.4) where κiare constants to be determined. No modification is introduced in equation (9.2). 10.2.2 Non-hysteretic term parameter identification The modification of the model implies a modification of the identification methodology. Similar to Section 9.1.3, two inputs x(t) and x(t) + qthat 106 10.2 The modified model and the corresponding identification methodology differ by a constant qare applied to the piezoelectric actuator. Then, the obtained asymptotic outputs ¯ Φ1(τ) and ¯ Φ2(τ) can be written: ¯ Φ1(τ) = N X i=1 κixi(τ) + κw¯w1(τ) (10.5) ¯ Φ2(τ) = N X i=1 κi(x(τ) + q)i+κw¯w2(τ) (10.6) where τ∈[0, T]. Note that we have ¯w1(τ) = ¯w2(τ),¯w(τ) from Theorem 2. Subtracting (10.6) from (10.5) it follows: ¯ Φ2(τ)−¯ Φ1(τ) = N X i=1 κih(x(τ) + q)i−xi(τ)i, N−1 X j=0 gjxi(τ),∀τ∈[0, T] (10.7) where gjare constant coefficients. Expanding the terms of equation (10.7) and rearranging we get:          g0 g1 g2 g3 . . . gN−1          =                        q q2q3q4. . . qN 0 2q3q24q3. . . N N−1qN−1 0 0 3q6q2. . . N N−2qN−2 0 0 0 4q . . . N N−3qN−3 . . .. . .. . .. . ..... . . 0 0 0 0 . . . N N−kqN−k . . .. . .. . .. . ..... . . 0 0 0 0 . . . N 1q                        ×          κ1 κ2 κ3 κ4 . . . κN          (10.8) which is equivalent to: 107 10. Adaptation of the Bouc-Wen model for the modeling and validation of a piezoelectric actuator Figure 10.6: Points used to determine the Bouc-Wen model parameters κw, n,ρand σ. In bold filtered experimental data. In grey fitted data for the computation of the derivatives. Table 10.5: Bouc-Wen model parameters Nn ρ κwσ 1 1.27 0.00893 5.08e-006 0.74 2 1.12 0.0047 9.22e-006 0.812 3 1.12 0.00463 9.35e-006 0.815 4 1.12 0.00461 9.37e-006 0.815 114 10.4 Conclusion model can be obtained from equation (10.15) as w(0) = y(0) −PN i=1 κivi(0) −y0 κw (10.18) where y(0), κi,v(0), y0and κware available. Figures 10.7(a) and 10.10(a) give the output of the model (10.15) for N= 1 and N= 2, along with the experimental output of the actuator. It can be observed that the model matches better the experimental data for N= 2. This conclusion can also be drawn from Figures 10.7(b) and 10.10(b), where the difference between the model and the experimental output is plotted for N= 1 and N= 2. It can be observed that after a transient phase, the error is smaller for N= 2. The same conclusions can be drawn from the displacement-voltage plot of Figure 10.9. Other experiments for N≥3 show that the behavior of the model for such values of Nis not significantly different from that of N= 2. 10.4 Conclusion The chapter has focused on the modeling of a piezoelectric actuator using a modified version of the hysteresis Bouc-Wen model. The modification consists in representing the non-hysteretic part of the model as a degree N polynomial instead of a linear relationship. The results for different values of N have been computed and compared with the real displacements of the actuator. The modified model has proven to match better the experimental data for N > 1. 115 10. Adaptation of the Bouc-Wen model for the modeling and validation of a piezoelectric actuator (a) Experimental output and model output for N= 1 and N= 2. (b) Model error for N= 1 and N= 2. Figure 10.7: Model response to a sinusoidal input. 116 10.4 Conclusion Figure 10.8: Excitation voltage Figure 10.9: Displacement - Voltage plot of the response to a sinusoidal input 117 10. Adaptation of the Bouc-Wen model for the modeling and validation of a piezoelectric actuator (a) Experimental output and model output for N= 1 and N= 2. (b) Model error for N= 1 and N= 2. Figure 10.10: Model response to a random input. 118 10.4 Conclusion Figure 10.11: Excitation voltage 119 10. Adaptation of the Bouc-Wen model for the modeling and validation of a piezoelectric actuator 120 Chapter 11 Control of a piezoelectric actuator considering the hysteresis This chapter deals with the modeling and control of a piezoelectric actuator. The main challenge for the control design is the presence of hysteresis. This nonlinearity is represented in this chapter using the Bouc-Wen model and a time-varying PID controller is designed for micropositionning purpose. The performance of the controller is tested using numerical simulations and experimentally. We consider the problem of micropositionning using a piezoelectric actuator. This problem has spurred much interest in the current literature. A robust controller is employed in [6] to control a piezoelectric bimorph actuator using the Bouc-Wen model. In [24] a piezoelectric actuator is modeled with neural networks and controlled with a variable structure control system. In [59], the controller uses information of the charge instead of the voltage for the control of position. This technique takes advantage of the reduced hysteresis between the displacement and the electrical charge, but presents some difficulty for the measurement of the charge. Since the piezoelectric device is represented in this work using the Bouc-Wen model, the results of [31] are used and improved for the control of the piezoelectric element. In 121 11. Control of a piezoelectric actuator considering the hysteresis [31], a second-order mechanical system that includes a Bouc-Wen hysteresis is considered for control purposes. The control objective is to guarantee the global boundedness of all the closed loop signals, and the regulation of both the displacement and the velocity of the device to zero. This objective is achieved using a simple PID controller. However, the main drawback of this controller is that the equilibrium point of the closed loop system is not robust vis-`a-vis perturbations which is undesirable in practice. The main contributions of this chapter are the following: •We present a new control law which is a time-varying PID that guarantees that the equilibrium point of the closed loop is robust to perturbations. •This control law is tested in numerical simulations and experimentally using a piezoelectric actuator. The main advantage of the proposed control law over other existing control schemes, is that it is simple to implement in an industrial context. 11.1 Background results. PID control of a Bouc-Wen hysteresis We consider the second order mechanical system described by: m¨x+c˙x+ Φ(x)(t) = u(t),(11.1) with initial conditions x(0), ˙x(0) and excited by a control input force u(t). The output restoring force Φ is assumed to be described by the normalized Bouc-Wen model [29]: Φ(x)(t) = κxx(t) + κww(t),(11.2) ˙w(t) = ρ˙x(t)−σ|˙x(t)||w(t)|n−1w(t)+(σ−1) ˙x(t)|w(t)|n(11.3) with an initial condition w(0). The parameters n≥1, ρ > 0, σ≥1 2, κx>0, κw>0, m > 0 and c≥0 are unknown. The range of the 122 11.1 Background results. PID control of a Bouc-Wen hysteresis parameters corresponds to the Class I Bouc-Wen model which is stable, asymptotically dissipative and thermodynamically consistent [29]. The displacement x(t) and velocity ˙x(t) are available through measurements, but the signal w(t) is not. Let yr(t) be a (known) smooth and bounded reference signal whose (known) smooth and bounded derivatives are such that limt→∞ yr(t) = limt→∞ ˙yr(t) = limt→∞ ¨yr(t) = limt→∞ y(3) r(t) = 0 exponentially. This means that there exist some constants a > 0 and b > 0 such that y(i) r(t)≤ae−bt for t≥0 and i= 0,1,2,3. The control objective is to globally asymptotically regulate the displacement x(t) and velocity ˙x(t) to the reference signals yr(t) and ˙yr(t) preserving the global boundedness of all the closed loop signals; that is x(t), ˙x(t), w(t) and u(t). We assume the following: Assumption 2 The unknown parameters lie in known intervals. That is we have m∈[mmin, mmax]with mmin >0,c∈[0, cmax],κx∈(0, κxmax ], κw∈(0, κwmax ],σ∈1 2, σmax,ρ∈(0, ρmax]. Note that the unknown structure parameter n≥1 is not required to lie in a known interval. The problem of controlling the system (11.1)-(11.3) has been treated in [31], where it is demonstrated that a PID control insures that the displacement and velocity errors tend to zero. Introduce the variables: x1(t) = x(t)−yr(t), x2(t) = ˙x(t)−˙yr(t), x0(t) = Zt 0 x1(τ)dτ (11.4) and choose as a control law the PID controller: u(t) = −k0x0(t)−k1x1(t)−k2x2(t) (11.5) where the ki’s are design parameters. Then we have [31]: 123