Modelling and control of the vibration of a rotatory disk

Master’s Thesis Double Master’s degree in Industrial Engineering and Automatic Control and Robotics Modelling and control of the vibration of a rotatory disk REPORT September 22, 2021 Author: Xavier Sánchez Botello Supervisors: Francesc Xavier Escaler Puigoriol Ramon Costa Castelló Delivery: Spring 2021 Escola Tècnica Superior d’Enginyeria Industrial de Barcelona Modelling and control of the vibration of a rotatory disk p. 1 Summary Complex systems submerged into water, such as hydraulic turbines, are subject to extreme and continuous forces that produce vibration during transient and off-design operations, which yield a decrease in efficiency, unit performance and lifetime. For that reason, as a part of a EU-funded AFC4Hydro project, the IFLUID’s research group from UPC has built a test rig in the laboratory with a rotating and submerged disk to investigate the dynamic behavior of that structure. This test rig is equipped with different sensors and actuators that permit to extract representative indicators of the system. The aim of this project is to develop a controller based on the dynamic response of a disk to reduce the vibrations produced by an external excitation, which can come from the machine itself or from any external factor. To identify the structure response, a piezoelectric patch has been used to excite the disk, and several strain gauges have been used to measure the strain produced by the vibration of the disk and also to determine the Frequency Response Function (FRF) of the system. Then, based on that FRF, the controller has been designed to determine the adequate control signal to feed the same piezoelectric patch and counteract the disk vibrations induced by the external excitations during operation. To design this controller, a numerical model of the dynamic response of the system has been developed in ANSYS and validated experimentally in the laboratory. The reason for its use, instead of a complete vibrational testing, has been justified based on the amount of time and resources that are saved. The simulation has permitted to extract the FRF of the system both in air and when it is submerged into water. Afterwards, it has been identified a plant model, in a form of a state space model, with the same FRF than the simulation. Once the plant model has been validated experimentally using LabVIEW, it has been used to design, develop and test a controller using MATLAB to reduce the vibrations induced in that simple rotating structure. All in all, this project encompasses the accurate simulation of the disk submerged into the tank of water, using a 2-way Fluid-Structure Interaction (FSI) analysis in ANSYS, the identification of a plant model in MATLAB from that complex simulation, and the experimental validation of it, using LabVIEW. Finally, the design of different control approaches has been done using MATLAB to end up reducing the vibrations of the structure. Modelling and control of the vibration of a rotatory disk p. 3 Acronyms AFC Active Flow Control APDL ANSYS Parametric Design Language CFD Computational Fluid Dynamics DOB Disturbance Observer FE Finite Element FEA Finite Element Analysis FFT Fast Fourier Transform FRD Frequency Response Data FRF Frequency Response Function FRR Frequency Reduction Ratio FSI Fluid-Structure Interaction LHP Left Half-Plane LQG Linear Quadratic Gaussian LQR Linear Quadratic Regulator NI National Instruments ROM Reduced Order Model SHM Structural Health Monitoring SIMO Single-Input Multiple-Output SISO Single-Input Single-Output TDMS Technical Data Management Streaming WEEE Waste from Electrical and Electronic Equipment WP Work Packages Modelling and control of the vibration of a rotatory disk p. 5 Glossary ANSYS is a Finite Element (FE) analysis software used to simulate computer models of structures, electronics, or machine componenets for analyzing strength, toughness, elasticity, temperature distribution, electromagnetism, fluid flow, and other attributes. ANSYS Fluent is a Computational Fluid Dynamics (CFD) software known for its advanced physics modeling and renowned for industry leading accuracy. ANSYS Mechanical is a FE solver for structural engineering with structural, thermal, acoustic, transient and nonlinear capabilities to improve the modeling in ANSYS. ANSYS Simplorer is an ANSYS multidomain simulation tool for complex systems, which utilizes models at various levels of detail, going from heavily simplified approximations to full high fidelity transient simulations. ANSYS Workbench is a new modern interface of ANSYS for solving engineering problems, in which you can import models from CAD systems, create the models using the Design- Modeler platform, and simulate the systems performing FE analysis. Bode diagram is a graph of the frequency response of a system. It is usually a combination of a Bode magnitude plot, expressing the magnitude of the frequency response, and a Bode phase plot, expressing the phase shift. LabVIEW stands for Laboratory Virtual Instrument Engineering Workbench, and is a systemdesign platform and development environment for a visual programming language from National Instruments (NI). Macro is a single instruction that expands automatically into a set of instructions to perform a particular task. MATLAB is an abbreviation of "matrix laboratory", is a proprietary multi-paradigm programming language and numeric computing enviornment developed by MathWorks. Simulink is a MATLAB-based graphical programming environment for modeling, simulating and analyzing multidomain dynamical systems. Its primary interface is a graphical block diagramming tool and a customizable set of block libraries. Skewness is defined as the difference between the shape of the cell and the shape of an equilateral cell of equivalent volume. Highly skewed cells can decrease accuracy and destabilize the solution. Step response is the time evolution of the outputs of a general system when its inputs change from 0 to 1 in a very short time. p. 6 Report Modelling and control of the vibration of a rotatory disk p. 7 Contents Acronyms 2 Glossary 3 1 Preface 15 2 Introduction 17 2.1 Objectives ......................................... 19 2.2 Scopeoftheproject.................................... 20 3 Experimental setup 21 4 Numerical modelling of the structure dynamic response 27 4.1 Dynamic response of the structure in air . . . . . . . . . . . . . . . . . . . . . . . . 27 4.2 Dynamic response of the structure in water . . . . . . . . . . . . . . . . . . . . . . 32 4.3 Comparison of the dynamic response of the structure in air and water . . . . . . 38 5 Extraction of the controller plant 41 5.1 Plant identification of the disk in air . . . . . . . . . . . . . . . . . . . . . . . . . . 42 5.1.1 Modal analysis of the disk in air . . . . . . . . . . . . . . . . . . . . . . . . 43 5.1.2 Harmonic analysis of the disk in air . . . . . . . . . . . . . . . . . . . . . . 46 5.1.3 Identification of the plant from the frequency response in air . . . . . . . . 51 5.2 Plant identification of the disk in water . . . . . . . . . . . . . . . . . . . . . . . . . 57 5.2.1 Modal acoustic analysis of the disk in water . . . . . . . . . . . . . . . . . 57 5.2.2 Acoustic harmonic analysis of the disk in water . . . . . . . . . . . . . . . 58 5.2.3 Identification of the plant from the frequency response in water . . . . . . 61 5.3 Comparison of the plant of the disk in air and water . . . . . . . . . . . . . . . . . 65 6 Experimental validation of the model 67 6.1 Validation of the disk in air . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 6.2 Validation of the disk in water . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 7 Design of the controller 77 7.1 Controllability and observability of the system . . . . . . . . . . . . . . . . . . . . 78 7.2 Control approach 1: LQR controller + observer . . . . . . . . . . . . . . . . . . . . 79 7.2.1 LQR controller on the plant model of the disk in air . . . . . . . . . . . . . 81 7.2.2 LQR controller on the plant model of the disk in water . . . . . . . . . . . 84 7.2.3 Implementation of the Kalman filter observer . . . . . . . . . . . . . . . . 85 7.2.4 LQG controller (LQR controller + Kalman filter observer) on the plant modelofthediskinair ............................. 87 7.2.5 LQG controller (LQR controller + Kalman filter observer) on the plant modelofthediskinwater............................ 89 7.2.6 Test the control approach 1 against external disturbances . . . . . . . . . . 91 7.3 Control approach 2: cancellation of poles and zeros . . . . . . . . . . . . . . . . . 92 7.3.1 Test the control approach 2 against external disturbances . . . . . . . . . . 99 7.4 Comparison between approaches . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 8 Disturbance rejection control 105 Modelling and control of the vibration of a rotatory disk p. 15 1 Preface Nowadays, there is an increasing trend to rely on renewable energy to produce electricity, aiming to combat climate change and to reduce the greenhouse gas emissions that come from the burning of fossil fuels. In order to reduce CO2emissions and local air pollution, which also affects human health, the world needs to rapidly shift towards low-carbon sources of energy such as renewable technologies. In fact, the growth of renewable energies aims to fully decarbonize the power sector by 2035, and achieve a net-zero carbon emissions by 2050. Renewable energy, often referred to as clean-energy mix, comes from natural renewable sources that are constantly replenished. Among all the renewable sources, the most popular and important ones are the solar, the wind, and the hydro energy. Solar and wind energy generation are very sensitive to the time and weather, as well as the season and the geographical location of the source. On the contrary, hydro-power plants are more robust to those changes because they can be created just building a dam or a barrier, and then exploiting the water level differences to generate electricity by releasing a flow of water from the reservoir, which spins a turbine that activates a generator to produce electricity. This energy is then used, among other applications, to regulate the electrical energy supplied when there is a peak in the required energy demanded and there is not enough solar or wind energy available. Therefore, it is of interest to increase the efficiency in extracting the hydraulic energy from the water and reduce the carbon footprint of the process. Hydraulic turbines, which are the core of these hydropower plants, are subject to extreme and continuous forces that, when controlled, improve unit performance, efficiency, and lifetime. In order to control this undesired forces, the EU-funded AFC4Hydro project is developing an Active Flow Control (AFC) system including an Structural Health Monitoring (SHM) that will improve operation in regimes outside the design window, by reducing pressure fluctuations, loads and vibrations. These improvements will lead to an expanded operation window, which will allow hydropower to contribute to system services of the electrical grid and to maintain a more natural flow in hydropower schemes, ultimately benefiting the economy, society and environment [10]. The AFC4Hydro project, which receives funding from the European Union’s Horizon 2020 research and innovation program, is divided in 6 different Work Packages (WP), each of them in charge of studying and analyzing a different part of the project. The "Universitat Politècnica de Catalunya" (UPC) is the Project Coordinator of the whole AFC4Hydro project. As part of the project, UPC has to develop, test and validate a SHM system at various scales and couple it to a controller developed in order to build the AFC system. To do that, the UPC has several research groups and departments involved in the project, which are: the "Centre d’Innovació Tecnològica en Convertidors Estàtics i Accionaments" (CITCEA-UPC), the "Institut de Robòtica i Informàtica industrial" (IRI) and the research group Barcelona Fluids & Energy Lab (IFLUIDS). This project represents a part of the work developed in UPC, that is centered on the study of the dynamic behavior of a rotating and submerged structure built in the laboratory by the IFLUIDS’s research group. As aforementioned, the aim of that test rig is to develop a SHM system to investigate the dynamic behavior of an hydraulic turbine, operating at off-design and at transient conditions. For that, different sensors have been installed in the structure to extract representative indicators of the structure vibration conditions. Modelling and control of the vibration of a rotatory disk p. 17 2 Introduction Structural vibrations are defined as continuous cyclic motions that can be measured and observed in an structure due to an excitation. When the structure is subjected to external forces, it tends to vibrate naturally at certain frequencies, known as the natural frequencies of the structure. These natural frequencies are dependent on the mass and stiffness of the structure, as well as its geometry. Due to that natural frequencies, when a structure is subject to a cyclic force whose frequency is equal or nearly equal to its natural frequency, it occurs a phenomenon called resonance. This resonance phenomenon makes the structure vibrate with larger amplitude than when the same cyclic force is applied at other frequencies, producing different motion patterns, which are known as structural mode shapes. The unwanted structural vibrations can cause fatigue or degrade the performance of the structure, decreasing their efficiency, producing cracks, heavy noise and in the worst case lead to even failure [11]. In the case of systems that transfer between different types of energy, such as the case of the hydraulic turbines, those vibrations can naturally emerge, so it is of interest to control that vibrations to avoid future damage to the structure. In order to control that structural vibrations, two different control techniques have been traditionally used: the passive control and the active control. Passive control methods are based on reducing the vibration by dissipating energy of the system without adding to it any external power. These type of methods include material damping enhancement, viscoelastic or frictional dampers, and various vibration absorbers and isolation schemes. The advantages of the passive approach are that the devices are relatively simple and cheap, and the system will always be stable since the control is realized by energy dissipation. However, since it produces fixed designs, such schemes might not be effective when the operating condition change. In addition, these schemes usually work well at the high-frequency region or within a narrow frequency range, but often have poor low-frequency performance. For that reason, active control methods are used to suppress the vibration of the structure, which are based on feedback/feedforward control schemes. A typical active control system consists of the plant, actuator(s), sensor(s), and the control electronics. The vibration control is then achieved by applying a secondary input to the structure, thereby modifying the system dynamic response to a desirable pattern. While active systems are generally more effective than passive methods, they have the disadvantages of being complicated and expensive, having the potential to destabilize the system, and being sensitive to system modeling error and uncertainties [12]. In the case study, active control methods are going to be used to control the vibration of the disk placed in the test rig mounted at the laboratory. Hence, the active control is going to be done using different piezoelectric patches mounted on the disk, acting as actuators, and strain gauges as sensors. As the output measured from the system will be the strain, the active control system to reduce the vibrations will also aim to reduce the strain produced by that vibrations, which is going to be monitored using strain gauges. In order to design that controller, it is needed to derive a dynamic model from the real test rig to know the behaviour of the disk. Analytically deriving the model of a complicated mechanical system such as the rotating disk is not a simple task. One way to obtain that model could be to perform a vibration analysis of the structure in order to predict its behavior under different dynamic conditions. However, since real systems are usually quite complicated when p. 18 Report viewed in detail, an exact physical characterization of that system is often impossible because it requires a large amount of sensors and the testing of many operating conditions. Thus, simplifying assumptions must be done to reduce the system to an idealized version whose behavior approximates that of the real system, process which is called modeling [13]. The modeling of the disk is done numerically simulating the structure and obtaining the different natural frequencies and a frequency response function of the structure, which is the structural response measured over a frequency range containing several resonance frequencies. The test rig also includes a tank of water to submerge the disk in water, so it has to be taken into account that the dynamic response of the disk will be affected by the presence of the surrounding water when the disk is submerged inside the tank of water. In this project, it is first presented the geometry and components of the test rig in Section 3, where the experimental setup of the system is described. On it, the different actuators and sensors used to excite and monitor the structure are detailed, together with the material they are made of and some useful properties that we need to introduce into the simulation. Then, the preliminary study of the project is focused on developing a numerical simulation of the structure to obtain its dynamic response when it is excited with a piezoelectric patch as an actuator, presented in Section 4. Concretely, the structure is modeled when it is rotating in air and when it is rotating inside the tank of water, where more complex factors have to be taken into account, such as the Navier-Stokes equation to describe the motion of the fluid. For that reason, a 2-way FSI simulation is needed to model the disk inside the tank of water, where the mechanical interaction of the disk with the surrounding water of the tank has to be modeled using both ANSYS Mechanical and ANSYS Fluent in a coupled scheme. Due to its complexity, this coupled simulation is detailed step by step in the Appendices of the project. After being able to simulate the transient dynamic response of the disk when excited with a voltage, a way to extract a plant model from the ANSYS simulation was studied and then presented in Section 5, based on identifying a plant model with a similar frequency response than the numerical one. This will determine the physical properties of the system that have to be taken into account in the plant model of the system, used to design the controller in both air and water scenarios. To verify that the plant model obtained from the simulation can be used to describe the real dynamics of the structure of the test rig, in Section 6 it is performed the validation of that model by programming in LabVIEW the different sensors and actuators that are placed on the disk. Then, the validation is performed comparing the frequency response obtained using the sensors in the laboratory with the one obtained using the numerical model. Finally, in Sections 7 and 8 the controller is designed to reduce the vibrations of the disk, which is based on reducing the amount of strain received by the sensors and computing the corresponding output voltage to the piezoelectric patch to reduce its effect. In those sections, different control approaches are presented, which allow to track a constant "safe" strain reference that ensures that the disk is not vibrating in the zone around the strain gauge. Additionally, in the latter section a disturbance rejection control based on a disturbance observer approach is designed to cancel any type of disturbances affecting the control action or the model of the system. Modelling and control of the vibration of a rotatory disk p. 19 Attached to the report are also provided different Appendices to explain in detail some concepts that are of interest for this project, such as the numerical modeling of a piezoelectric body, presented in Appendix A, or the material and geometry of the numerical ANSYS model, presented in Appendix B. In Appendices C and D the steps performed to simulate the system using ANSYS Mechanical and ANSYS Fluent are detailed, whereas the coupling of both systems to perform the 2-way FSI simulation is detailed in Appendix E. Then, in Appendix F is detailed how to introduce different non-standard voltages into the simulated piezoelectric element. Appendix G is devoted to explain the procedure to extract the plant model from the ANSYS simulation and import it into MATLAB, together with the code and other relevant information of this process. Finally, in the last Appendix H the MATLAB code used to program the controllers is presented, together with the suitable comments to understand each line of code. 2.1 Objectives The main objective of the project is to design a controller to reduce the vibrations of a submerged rotating structure. To do it, first we need to obtain a numerical model of the system that represents its dynamic behavior under the same operating conditions as in the laboratory. Then, this model has to be validated using experimental data to ensure that it represents the dynamic behavior of the test rig. For that, piezoelectric actuators will be used to excite the structure and strain gauges will be used to measure its response. Finally, this numerical model will be used to identify a plant model and design its corresponding controller to reduce the vibrations and strain of the disk. Therefore, the partial objectives that are needed to achieve the main objective of the project are the following ones: •Define all the relevant physical properties of the real structure needed to build a simplified numerical model in ANSYS. •Integrate the piezoelectric actuation model in the ANSYS simulation. •Obtain numerically the dynamic response of the system in air and water. •Determine numerically the natural frequencies and the corresponding mode shapes of the system when it is in air and under water. •Obtain numerically the frequency response of the structure in air and water. •Identify a plant model of the system based on the frequency response data obtained from the numerical simulations. •Validate the numerical plant model using the sensors and actuators of the laboratory and the real test rig. •Apply different control techniques to design a controller to minimize the vibrations of the system. •Validate the effectiveness of the controller. p. 20 Report 2.2 Scope of the project In this project different fields will be covered, going from the modelling of a test rig built in the laboratory to the experimental testing. For the modelling it will be investigated the best way to obtain and develop a numerical model for performing a simulation, whereas in the experimental testing we will program the sensors installed in the disk and analyze its results. Finally, different controllers will be designed using the suitable control approach to achieve the control specifications. The construction and validation of the test rig has already been done by the IFLUIDS research group, it is only measured and numerically simulated what is already built. Moreover, the sensors used to acquire experimentally the data and the piezoelectric actuators are provided by them, so they are not chosen specifically for this project. Modelling and control of the vibration of a rotatory disk p. 21 3 Experimental setup The system considered in this project is a rotating structure constructed in the laboratory of UPC, which is basically a thin disk attached to a long shaft that is placed on a test rig that holds it and allows it to rotate. The test rig, as it can be seen in Figure 3.1, is composed with different metallic beams that can be easily adjusted at different positions by unscrewing and screwing some bolts. This metallic beams are used to support the bearings of the central shaft of the disk, as well as the rotating motor of the system. Then, once the beams are placed on the suitable position to support the bearings of the shaft, the disk can be placed in a specific depth inside the tank of water by means of softening and tightening the bearings supporting the shaft. Moreover, the test rig is also equipped with a DC motor coupled to the shaft of the disk using a belt, which allows to rotate the disk at different speeds. Figure 3.1: Experimental setup of the test rig mounted in the laboratory. The tank of water is a cylinder of plexiglass with an internal diameter of 488 mm and a height of 600 mm, which is transparent to allow the user to see the interaction between the disk and the water. The rotating disk is a thin circular plate of 420 mm of diameter and a thickness of 2 mm. The through hole of the disk, which also passes through the shaft, has a diameter of 30 mm and a thickness of 6 mm. The shaft that supports the disk has a length of 705 mm, and it has two bearings of 35 mm of height each one, that enables the support of the shaft to the structure of the test rig. Both bearings are separated 168 mm from each other using the metallic beams of the test rig. In Figure 3.2 it is shown a drawing scheme of the structure made in SolidWorks, where the relevant dimensions (in mm) are shown. p. 22 Report Figure 3.2: SolidWorks drawing of the disk submerged 150 mm inside the tank of water. Dimensions in mm. In the particular case of the figure, the disk is submerged 150 mm inside the tank of water, so depending on the depth of the disk inside the tank of water the relative position of the bearings will change accordingly, as well as the depth of the disk inside the tank. The coupling between the shaft and the disk is done using a set of bolts that are directly screwed to the shaft from the bottom of the disk. Therefore, instead of considering the complex geometry made of several bolts used in the real structure, the geometry is simplified considering that all the threads of the bolts form a cylinder of 63.60 mm of diameter and 11 mm of height. The test rig is also equipped with systems to generate mechanical and fluid-dynamic excitations, together with measuring systems to check the deformation and vibration of the structure. As it is needed to collect the data from the sensors in real time, the disk and the shaft have a through hole with a slip ring on the top of the shaft to allow the transmission of power and electrical signals from the stationary data receiver device to the rotating structure of the disk. The additional elements used in the disk are piezoelectric transducers, working as actuators, and different sensors such as strain gauges, lasers, accelerometers, pressure sensors, and fiber optics. The piezoelectric transducers used in the system are bonded to the disk surface using a water resistant epoxy adhesive. The piezoelectric transducers used in the test rig are the P-876.A12 and P-876.SP1 DuraAct patches commercialized by PI Ceramic, both of them consisting of a PIC255 active layer sandwiched between two soft thin encapsulating polymer materials, together with Modelling and control of the vibration of a rotatory disk p. 23 two electrodes, as presented in Figure 3.3a. This polymer coating serves as electrical insulation for the active layer and is used to apply prestress to the active layer in order to avoid cracks, which increases the allowable curvature radius of the PZT ceramics. The overall dimensions of the P-876.A12 transducer (including the encapsulating material) are 61 mm x 35 mm x 0.5 mm, while the encapsulated active piezoceramic layer is 50 mm x 30 mm x 0.2 mm. In Figure 3.3b it is shown the drawing of the P-876.A12 patch with the dimensions in mm, together with the tolerances provided by the manufacturer [1]. (a) Layers of the piezoceramic element. (b) Dimensions of the piezo, in mm. Figure 3.3: DuraAct transducer P-876.A12 design, geometry and dimensions. Images from [1]. As it is not possible to know exactly the thickness of the epoxy glue used to bond the piezo to the disk, it is approximated to a value of 0.2 mm. Having this in consideration, we end up with a section of 5 layers of the piezoceramic element, presented in Figure 3.4. The outer polyimide of the patch is made of Kapton, which is a suitable coating material due to its good bonding property and electrical insulation, allowing the patch to be submerged in water. This Kapton layer has a total thickness of 0.3 mm and it is surrounding the active layer, so 0.15 mm of Kapton it is considered above and below the active layer, which has a total thickness of 0.2 mm. Considering the approximated value of the glue of 0.2 mm between the patch and the disk, the geometry ends up with a soft layer of 0.35 mm between the active layer and the disk, used as a medium to represent the combined contributions of the glue and the encapsulating Kapton shell [9]. Figure 3.4: Schematic of the section of a P-876.A12 patch. To ascertain the amount of deformation caused by an external force applied to the test rig, some strain gauges are used to measure the strain, which is a dimensionless measure (m/m) that relates the elongation of the material when applied a force ∆L, with the original length L, obtaining that the strain is computed as ∆L/L. The aforementioned strain gauges are the sensors that change its electric resistance proportionally to the deformation of the measurement object, so they measure strain by means of measuring the resistance change. p. 30 Report The stress obtained, which has a maximum value of 2.4054 MPa just below the piezoelectric patch, strongly depends on the type of piezoelectric used. For this simulation, the P-876.A12 model of the DuraAct patches from PI ceramics is used, which has an operating voltage range that goes from -100 to 400 V, having a blocking force of 265 N. Thus, this stress value will change accordingly depending on the dimensions of the patch used in the laboratory and its corresponding blocking force, which is defined as the force that the actuator exerts when it is blocked from moving. Figure 4.3: Equivalent Stress of the disk in air after applying 100 V in the piezoelectric patch. When this stress induced by the patch is applied to the disk, it produces a change in the shape of the object, which is quantified using the strain measure (m/m). This measure can be monitored in the laboratory using the different strain gauges of 1 and 3 axes available, so it is of interest to know the strain distribution of the disk when it is rotating, at it will give an idea of the zones where the strain is maximum and it could be optimum to place a strain gauge to see a significant change of the strain when different voltages are applied to the patch. In Figure 4.4 it is shown this strain distribution, which is maximum around the piezoelectric patch (in the bottom face of the disk) and it is also significantly big in the top face of the disk, aligned with the position of the patch. Figure 4.4: Equivalent Elastic Strain of the disk in air after applying 100 V in the piezoelectric patch. Modelling and control of the vibration of a rotatory disk p. 31 Both stress and strain distributions are very similar, as they have a linear relation between them represented by the Hooke’s Law, that relates both stress and strain measures with the Young’s modulus of elasticity of the structure, denoted as E. This modulus of elasticity is essentially a measure of the stiffness of each material, and is one of the factors used by ANSYS to calculate a material’s deflection under load. As all the different material properties are introduced in the simulation, ANSYS is also able to compute this deflection of the structure under the load applied by the piezoelectric, as presented in Figure 4.5, where the total deformation of the structure is presented (in m). Figure 4.5: Total deformation of the disk in air after applying 100 V in the piezoelectric patch. More than the total deformation of the structure, it is of more interest its vertical displacement when the voltage is applied, as it will help to understand how the structure will change its shape under a constant voltage applied to the patch. In Figure 4.6 the directional deformation in the Z axis is shown, where very small displacements are present, as expected. Figure 4.6: Directional deformation in the Z axis of the disk in air after applying 100 V in the piezoelectric patch. Looking at the vertical displacements it can be deduced the bending direction of the patch, together with the slightly deformation of the disk, whose tip nearer to the patch moves upwards with a total maximum displacement of 13.53 µm when the voltage is applied. On the other hand, the diametrically opposed tip to the patch moves downwards with a negative Z displacement of -2.3963 µm. The shaft, even though it has a bluish colour, has restricted its vertical displacement and acts as a fixed point, so it has vertical displacement of 0 m. p. 32 Report 4.2 Dynamic response of the structure in water The test rig is designed to be able to submerge the disk inside the transparent Plexiglas tank of water and see its hydraulic performance when rotating it at different speeds. Therefore, a dynamic response of the structure when it is rotating inside the water is also of interest. Different from the in-air simulation, the water induces a certain level of complexity into the problem, as it has to be taken into account that the pressure distribution of the fluid will cause an extra mechanical stress at each time instant, and it will reduce the deflection of the disk due to the damping effect of the water. The system considered in this case is the scenario presented in Figure 4.7, whose numerical model includes now the tank of water, as it is shown in the Figure 4.7b. (a) Disk submerged in water. (b) Numerical modelling of the disk under the water. Figure 4.7: Geometry considered for simulating the disk in water. One way to deal with this water around the disk is to apply the pressure obtained from a CFD simulation as a boundary condition in the previous mechanical FEA simulation of the structure, which is called a 1-way coupled FSI, since the deflection of the solid is not fed back into CFD. On the other hand, if we want to capture the influence of the deformed solid on the hydraulic performance, the deformation of the disk has to be feed back into the CFD solution to close the loop and obtain an improved solution closer to the reality. This second approach is called the 2-way coupled FSI, and it obtains a more reliable model with an environment closer to the laboratory conditions [16]. Therefore, this approach is the one used in order to obtain the numerical model of the structure submerged in water and simulate its dynamic response when the piezoelectric element is excited with a certain voltage. The steps that are needed to perform this 2-way coupled FSI analysis are the following ones: first, a Transient Structural FEA in ANSYS Mechanical is needed to calculate the mechanical stresses and deflections of the disk. Then, a Fluid Flow CFD study in ANSYS Fluent is needed to solve the pressure-based Navier-Stokes equations and obtain the pressure distribution of the water around the disk. Finally, once each system is modeled in the corresponding simulation, Modelling and control of the vibration of a rotatory disk p. 33 the 2-way FSI analysis has to couple both models to coordinate and synchronize the solution execution. The general workflow of the system is presented in the Project Schematic of Figure 4.8, where data flows downstream (from top-to-bottom and from left-to-right across systems). The square connectors between two cells indicate that the port in both of them are shared together, while the round connector shows that the results from the source cell are transferred to the target cell. Figure 4.8: Project Schematic of the 2-way FSI used in ANSYS Workbench. In the project schematic it can be seen that a common "Engineering Data" and "Geometry" systems are used for both Transient Structural and Fluent simulations to define the material properties and the geometry of the model to simulate, which are detailed in Appendix B. The geometry model considers that the disk is submerged 150 mm under the water, whose free surface is separated 47 mm from the bearings that support the shaft. The solid used to model the reservoir of water is a cylinder representing all the volume occupied by the water, so it is needed to subtract all the structural solids from it. The next step is to define the Transient Structural simulation, whose configuration relating the piezoelectric body and the physical boundary conditions is analogous to the presented in the simulation in air, detailed in Appendix C. However, in this case it has to be specified which part of the solid is in contact with the water, as will be the responsible to send the displacement obtained in the mechanical analysis to the fluid analysis, and then introduce into the mechanical simulation the force data from that fluid analysis. This additional boundary condition is the Fluid-Solid Interface, explained in the subsection C.6. The final step is to define the environment of the incompressible fluid flow in ANSYS Fluent, which has the capacity to study the dynamic response of the water by solving the conservation equations of mass and momentum. As it is a complex simulation with a lot of parameters to be taken into account for a correct modelling, all the configuration used in the Fluent environment is detailed step by step in the Appendix D. At the end, both simulations are coupled using the "System Coupling" block, responsible to handle the data transmission at each time step between the systems. In Figure 4.9 it is shown how the different data is being transferred between the simulations at each coupling iteration, as well as the convergence chart of this transfer. Again, the detailed configuration of this data transfer in the "System Coupling" block is presented in Appendix E. p. 34 Report Figure 4.9: Outline of the System coupling block used in the 2-way FSI. As it has been defined a rotational velocity of 6.28 rad/s, a total simulation time of 1 second is fixed to obtain the dynamic response of the disk when it rotates one turn. Different voltages have also been tested on the piezoelectric patch to see the effect that it has on the results. Concretely, a fixed voltage value, a step and a sinusoidal voltage has been tested, which are introduced in the simulation using the APDL commands presented in Appendix F. However, a constant voltage of 100 V has been used as the "Voltage" boundary condition of the final simulation, as it is introduced directly with the PiezoAndMems extension. By solving the coupled system, it can be obtained the resulting mechanical forces applied to the disk when it is submerged into the water and rotating at a specific speeds, considering the constant voltage applied to the piezoelectric patch. On the other hand, in the Fluent simulation it is obtained the hydraulic behavior of the water of the tank under the influence of the rotating disk. The deformation of the disk at the time instant of t= 1 s, when it has turned one round, is presented in Figure 4.10. As expected, the zone of the disk with higher deformation is the extreme of the disk that is closer to the patch, which is constantly exerting a strain when the voltage is applied. On the contrary, the center of the disk, which is attached to the shaft, is presenting almost a null deformation, as its movement is constrained with the bearings boundary condition. Modelling and control of the vibration of a rotatory disk p. 35 Figure 4.10: Total deformation of the structure submerged in water after applying 100 V in the piezoelectric patch. The amount of the previous movement that corresponds to the vertical displacement in Z axis is shown in Figure 4.11, which has a similar pattern as the total displacement. That means that the majority of the displacement produced in the disk is due to the force applied by the patch in the vertical direction. Figure 4.11: Directional deformation in the Z axis of the structure submerged in water after applying 100 V in the piezoelectric patch. In the test rig, different strain gauges are used to detect and localize damage in the SHM system while the disk is rotating inside the water, as they are water-resistant and very sensitive to detect local damage. However, it is needed to know the optimal location to place those strain gauges in order to measure the maximum amount of strain. Hence, the 2-way FSI simulation can be used to know the strain distribution of the whole disk at each time instant and decide the optimal location to place those strain gauges. The strain distribution is obtained at each time-step of the Transient Structural results, obtaining the distribution shown in Figure 4.12, in which both faces of the disk are shown. We can see that the strain is maximum in the bottom face of the disk, where it is placed the piezoelectric patch, as it is constantly applying a stress to the disk. On the other hand, in the other side of the disk, p. 36 Report in the top face, the strain is big enough to see the influence of the piezoelectric patch, but lower than the other face because there is not a direct contact between the patch and the disk. Then, as long as we go move away from the location of the patch, the strain decreases to a negligible value. Figure 4.12: Equivalent Elastic Strain of the structure submerged in water after applying 100 V in the piezoelectric patch. As well as in the modelling of the structure in air, it is also obtained the equivalent stress of the structure when submerged into water, presented in Figure 4.13. On it, we can also see that there is a strong relation between the strain and the stress distribution of the disk, which are linearly related with the aforementioned Young’s modulus. Figure 4.13: Equivalent Stress of the structure submerged in water after applying 100 V in the piezoelectric patch. Regarding the Fluent simulation, the results are shown using the CFD-Post application, which allows multiple solution datasets to be loaded simultaneously. This application is used to merge the results obtained from the Transient Structural and the Fluent simulation, being able to visualize the synchronized 2-way FSI results in a compact manner. Modelling and control of the vibration of a rotatory disk p. 37 In Figure 4.14 it is shown the streamline distribution inside the tank of water when the disk is rotating, which clearly demonstrates that the particles of water inside the tank are rotating accordingly to the rotation of the structural simulated disk. This result matches with the approach explained in Appendices C and D of considering, on one side, the fictitious rotation in the Transient Structural simulation to account for the equivalent forces due to the rotation and, on the other side, the rotation equivalent in the walls that are in contact with the fluid, in the Fluent simulation. Even though this result seems to be very trivial, different convergence problems came up in the process of developing a 2-way FSI simulation with this rotation, as the incremental displacement sent at each iteration to the Fluent was not enough to simulate the rotation of the solid in Fluent. Moreover, this displacement yielded different convergence problems due to the mesh rotation, producing several errors when updating the mesh of the model. Figure 4.14: Streamlines of the tank of water after applying a constant voltage to the patch. Moreover, it can be seen that whenever the streamline that starts from the top face of the tank reaches the disk, it gets attached to it, and it starts rotating at the same speed of the disk. This linear velocity of the particles attached to the disk can be computed at each radius Rfrom the center of the disk using the equation v=ω·R, where ω= 6.28 rad/s. As an example, it can be seen that the maximum value of the velocity of the disk is 1.313 m/s, which corresponds to the linear velocity of any particle attached to the tip of the disk, at a radius of 210 mm. In order to see the effect that the rotating structure has on the pressure distribution inside the tank of water, a central plane has been created in the middle of the tank, as it is shown in Figure 4.15, in which it is depicted the pressure distribution of the water after 1s of simulation. It can be seen that in the bottom face of the disk there is a negative relative pressure, which makes a suction in that part of the tank of water. On the contrary, in the regions around the tip of p. 38 Report the disk, a higher pressure is applied due to the rotation of the disk, which makes the particles of water move due to this gradient of pressure inside the tank of water. However, this gradient of pressure inside the tank of water represents an insignificant change compared with the atmospheric pressure, which is 101,325 Pa. Figure 4.15: Pressure distribution of the central plane of the tank of water after applying a constant voltage to the piezoelectric patch. 4.3 Comparison of the dynamic response of the structure in air and water Once the two different scenarios have been studied separately, it has been obtained that the strain and stress distribution on the disk surface is very similar, as well as the deformed shape of the disk due to the voltage. Comparing in detail the results obtained from both scenarios, we can see that when the disk is submerged under the water the force received by the structure slightly increases due to its hydraulic pressure. This pressure is equivalent to a bulk stress, which is the kind of stress that represents a force pressing on every part of the submerged structure, applied from all directions. This force is perpendicular to the submerged surface an it slightly decreases the volume of the submerged object by a certain amount. Hence, this reduction of volume relative to the original one, produces an additional strain to the disk called bulk strain. This produces that the strain measured in the disk when it is submerged in the tank of water is bigger than in air, due to this extra factor considered [17]. However, as the differences between the maximum values of stress and strain in the disk are insignificant, it is concluded that the bulk stress effect produced by the water is almost negligible. Modelling and control of the vibration of a rotatory disk p. 39 Regarding the total deformation of the disk, it is bigger in the case of the disk being inside the tank of water, whereas the vertical displacement in the Z direction is almost equal in both cases. Theoretically, the vertical displacement of the disk inside the tank of water should be damped by the fluid, but in this case it is not modified because the mesh in the Fluent analysis is bigger than the displacements produced in the simulation by the patch. This produces that, numerically, the mesh of the fluid is not able to adapt to the small deformations of the disk, so it considers that the disk is barely deformed in the Z direction and only the rotation of the disk has a significant factor when it is submerged in water. Therefore, the increase of the total deformation of the disk when it is submerged in water is produced by the shear force that appears when the disk is rotating inside the water, rather than by the variation of the deformation in the Z direction. This shear force produces an additional tangential and axial deformation to the disk that increases the deformation received in the X and Y directions of the disk. To compare the effect of the shear stress to the deformation of the disk in X and Y directions, in Figure 4.16 is presented the different deformations in X direction obtained when the disk is in air and when it is in water. In this case, it can be clearly seen that in the case of having the disk rotating in air, this shear force does not affect to it, as the force produced by the patch is bigger than the friction of the disk with the air. On the other hand, when the disk is rotating inside the tank of water, it appears an extra deformation due to the friction of the disk with the surrounding water, increasing the deformation accordingly to the rotation. (a) Directional deformation in X direction in air. (b) Directional deformation in X direction in water. Figure 4.16: Comparison of the directional deformation in X direction. In the case of the deformation in the Y direction, the behavior is exactly the same as in the X direction, as presented in Figure 4.17. In this case, it can also be clearly seen that, when the disk is submerged in water, the deformation in the Y direction increases due to the shear force produced by the viscosity of the simulated surrounding water of the disk. p. 46 Report 5.1.2 Harmonic analysis of the disk in air To determine the response of the structure under a steady-state sinusoidal harmonic voltage excited at a given frequency, an harmonic analysis is performed. This harmonic analysis will give an idea of the strain amplitude measured at each natural frequency obtained from the modal analysis and will also be used as the model to identify the controller plant of the system. The most relevant assumptions made in an ANSYS harmonic analysis are that the entire structure has a constant or frequency-dependent stiffness, damping, and mass effects, and that all loads and displacements vary sinusoidally at the same known frequency (but they can be in a different phase). Governing equations of Harmonic Analysis in ANSYS The general equation of motion used by ANSYS to solve the harmonic analysis is the following one: [M]{¨u}+ [C]{˙u}+ [K]{u}={F}(5.5) where {u},{˙u},{¨u}are the vectors of nodal displacement, velocity and acceleration of the n discretized points in the time domain, respectively. [M]is the structural mass matrix, [C]is the structural damping matrix of the system, [K]is the stiffness matrix, and {F}is the force vector that expresses the force applied on each point in the time domain, which varies sinusoidally at a given frequency Ω[22]. All the points in the structure are assumed to be oscillating at the same frequency, so the complex force and displacement vector can be defined as the sinusoidal signal of equation (5.6), as a function of time tand a circular frequency Ω: {F}=FmaxeiψeiΩt=Fmax(cos ψ+isin ψ)eiΩt= ({F1}+i{F2})eiΩt {u}=umaxejφeΩt=umax(cos φ+isin φ)eiΩt= ({u1}+i{u2})eiΩt(5.6) where Ωis the excitation frequency at which the loading occurs, ψis the force phase shift that may be present if different loads are excited at different phases and φis the displacement phase shift that may exist if damping or a force phase shift is present. By substituting the expressions of equation (5.6) into the general equation of motion of (5.11), and differentiating the displacements to get velocity and acceleration, it can be obtained the harmonic equation of motion: (−Ω2[M] + iΩ[C]+[K])({u1}+i{u2})=({F1}+i{F2})(5.7) where {u1}={umax cos φ}is the real displacement vector, {u2}={umax sin φ}is the imaginary displacement vector, {F1}={Fmax cos ψ}is the real force vector and {F2}={Fmax cos ψ}is the imaginary force vector [23]. Modelling and control of the vibration of a rotatory disk p. 47 ANSYS has two ways of solving the harmonic equation of motion (5.7), depending on the solution method selected, which can be either the Mode Superposition Method or the Full Method. The Mode Superposition Method solves the harmonic equation in modal coordinates, as shown in Equation (5.8), where all the displacements are uncoupled and can be expressed as a linear combination of the mode shapes φj. ¨yj+ 2ωjξj˙yj+ω2 jyj=fj(5.8) Therefore, it is first needed to perform a modal analysis to determine these mode shapes φj and their corresponding natural frequency ωj. For each mode jit can also be extracted the modal coordinate yj(together with its acceleration ¨yjand velocity ˙yj), the fraction of critical damping ξjand the force applied in modal coordinates fj. However, the solution obtained is approximate, and it does not fully support nonzero displacements. On the other hand, the Full Method solves the matrix equation (5.7) directly in nodal coordinates, which is analogously to a linear static analysis, except that complex numbers are used. Although this method is more computationally expensive, it computes an exact solution and supports all types of loads including nonzero displacements, including the possibility to define the voltage on each face of the piezoelectric element. Hence, a Full Method is used for solving the harmonic analysis, which does not rely on mode shapes information and natural frequency. Simulation of an Harmonic Analysis in ANSYS To perform the harmonic analysis in ANSYS the block "Harmonic Analysis" is used, with the same geometry and material definitions than in the Modal analysis. In this case, the displacement of the bearings is also fixed, and the piezoelectric element body is defined using the PiezoAndMems extension. In order to be able to use voltage as an input of the harmonic analysis, the aforementioned Full Solution Method has to be selected as the "Solution Method" of the analysis. By doing so, we are able to define the desired harmonic input voltage into the POSITIVE face of the active layer and the 0 voltage into the NEGATIVE face. As well as in the modal analysis type, a constant structural damping coefficient of 0.02 is set, which is a typical value for metals. Solving the harmonic analysis with an input voltage of 100 V to the piezoelectric patch, the strain frequency response of each strain gauge can be obtained, as presented for the "Next" strain gauge in Figure 5.4. The solution of that analysis contain a graphical representation of the frequency response and the corresponding tabular data of the strain amplitude (m/m) and phase angle (◦) associated to each frequency (Hz). p. 48 Report Figure 5.4: Solution of the harmonic analysis for the "Next" strain gauge obtained in ANSYS, with an input voltage of 100 V. To be able to obtain a more accurate model around the natural frequencies, different sampling rates are used along the frequency range of 20-326 Hz, using in the most important zones a sampling rate of 0.001 Hz, and 0.2 Hz in the transition zones. This refinement of the frequency response has been done repeating the ANSYS simulation with more points around each natural frequency, to obtain a better resolution. In Appendix G.2 the importance of the refinement process is presented, where it is shown that a lot of information is lost with a poor quality frequency response signal. Moreover, in the Appendix it is shown the difference on the Frequency Response Data (FRD) obtained when using a structural damping value of 2% and when using no damping, producing that the peaks and valleys tend to a very high and unreal value. Once obtained a reliable frequency response for the disk in air, it is needed to export the tabular data of the frequency response into an excel file to further import it into MATLAB, where all the System Identification Toolbox functions are available to identify a model. An example of this post-processing of the data is explained in detail in Appendix G.1, together with the MATLAB Script 5 that it used to extract the frequency, amplitude and phase of the frequency response from the excel file, and then store it as a FRD in MATLAB. As explained in the Appendix G.1, this FRD vector can be stored in MATLAB using different commands depending on the toolbox used. Hence, using the System Identification Toolbox, the data is stored using an idfrd object, where the FRD of both strain gauges are encapsulated in the same model with 1 input and 2 outputs. As an additional information, the names and the units of the input and the outputs are defined, being able to differentiate between the behavior of the "Top" and "Next" strain gauges. Finally, to ensure that the MATLAB idfrd object has the same frequency response as the source data of ANSYS, it can be computed the Bode diagram of each output, presented in Figure 5.5. Modelling and control of the vibration of a rotatory disk p. 49 Figure 5.5: Strain frequency response imported to MATLAB, with an input voltage of 100 V. These strain frequency responses (m/m) have to be normalized considering the input voltage of the system (V) in the frequency domain, obtaining the Frequency Response Function (FRF). The FRF is defined as the ratio of the complex output amplitude to the complex input amplitude for a steady-state sinusoidal input, as it is shown in equation (5.9). In Appendix G.3 is studied the variability that may appear on the computation of this FRF when simulating the system using different input voltages. On it, it is presented the different strain frequency responses obtained using 4 different input voltages (1 V, 10 V, 50 V and 100 V), ending up demonstrating that using simulated data the resulting FRF is independent on the voltage applied to the piezoelectric patch. FRF =Magnitude ·ej·Phase =AmplitudeStrain ·ej·PhaseStrain AmplitudeVoltage ·ej·PhaseVoltage (5.9) Knowing that the amplitude of the input voltage in this particular case is 100 V, and that the phase of the voltage is 0◦at each frequency, the strain frequency responses can be transformed to obtain the FRFs presented in the following Bode diagram of Figure 5.6, where the natural frequencies obtained in the previous modal analysis are marked using a green dashed line. Theoretically, a peak on the magnitude of the FRF should appear at each natural frequency, as the system starts vibrating and the strain of some zones of the disk increases. However, depending on the mode shape it may not produce strain in the direction of the gauge, as happens with the torsional mode of f0 4= 116.14 Hz and the shaft bending mode of frequencies f0 9= 320.42 Hz. In those cases, the effect of the natural frequencies cannot be determined using strain gauges, but in the rest of the mode shapes the peak of the magnitude effectively matches with the natural frequency of the system. p. 50 Report Figure 5.6: Bode diagram of both strain gauges obtained using an input of 100 V. The peaks obtained in the strain FRF, using the harmonic analysis, are compared with the natural frequencies of the modal analysis in Table 5.1. It can be seen that the natural frequencies of both gauges obtained from the harmonic analysis are exactly the same, as both measure the vibration of the same structure. Moreover, comparing the values from both modal and harmonic analysis, we can see that they are almost the same, as they correspond to the solution of the harmonic equation of motion used by ANSYS. Table 5.1: Comparison modal analysis and harmonic analysis in air. Mode shape iModal fi(Hz) Harmonic "Next" fi(Hz) Harmonic "Top" fi(Hz) Error "Next" % Error "Top" % Mode shape 1 (1ND) 43.157 43.167 43.167 0.023 % 0.023 % Mode shape 2 (0NC) 49.947 49.952 49.952 0.010 % 0.010 % Mode shape 3 (2ND) 63.842 64.797 63.797 0.070 % 0.070 % Mode shape 4 (Torsional) 116.14 - - - - Mode shape 5 (3ND) 133.27 133.056 133.056 0.161 % 0.161 % Mode shape 6 (Shaft bending) 156.42 156.456 156.456 0.023 % 0.023 % Mode shape 7 (4ND) 232.93 232.9 233.1 0.013 % 0.073 % Mode shape 8 (1NC) 305.74 305.8 305.8 0.020 % 0.020 % Mode shape 9 (Shaft bending 2) 320.42 - - - - Mode shape 10 (1NC-1ND) 330.80 330.606 330.606 0.059 % 0.059 % Modelling and control of the vibration of a rotatory disk p. 51 It is important to mention that the shaft bending mode of frequency f0 6it is theoretically not possible to be obtained from the gauges because it does not produce any strain to the disk (as happens with the torsional mode of f0 4and with the second shaft bending mode of f0 9). However, it appears in the current harmonic analysis because it is, in fact, not a purely bending mode but a combination of a shaft bending mode with another mode that deforms the disk, as it can be seen in the previous Figure 5.3. The dominant deformation is induced by the bending of the shaft, but it is also slightly deforming the disk, which is why the strain gauges are able to notice that strain produced by the disk and it can be obtained a natural frequency from them. 5.1.3 Identification of the plant from the frequency response in air At this point a frequency response of the simulated system has been obtained, in which the input is the voltage applied to the piezoelectric patch, and the output is the strain obtained in the 2 strain gauges ("Next" and "Top"). Therefore, a model that presents a similar frequency response is going to be identified using MATLAB. The approach done is to try to adjust the biggest model as possible, in order to capture the full FRF obtained, and then perform a model order reduction to keep only the most dominant poles and zeros of the model. Analyzing in detail the Bode diagram of the Figure 5.6, it can be deduced that the poles of the system correspond to the peaks of the Magnitude plot, while the zeros of the system are the valleys of that plot. The peaks of both outputs, corresponding to the poles, are all aligned at the same frequency because the poles are intrinsic properties of the system and they do not depend neither on the output sensors yneither on the input actuators u. On the contrary, the zeros of a system depend on how inputs and outputs of a system are coupled to the states, so zeros can be changed by moving sensors and actuators or by introducing new sensors and actuators. Therefore, the valleys are not in the same position for the different output strains. The MATLAB System Identification Toolbox has several functions for constructing mathematical models of dynamic systems from measured input-output data. Using a linear model is sufficient to completely capture the system dynamics, so parametric models such as transfer function models and state-space models are used to identify the model. Some examples of this functions are presented below, where nprepresents the number of poles of the identified system, nzthe number of zeros and nkthe discrete delay of an input disturbance: •Identify a Transfer Function model: tfest(FRD,np,nz), where np≥nzto be causal •Identify Output-Error polynomial model: oe(FRD,[np,nz]), which uses white noise as an additive output disturbance •Identify a State Space model using ODE: ssest(FRD,nx), where nxis the number of states In addition, for comparing the different models, some fitting parameters are provided by the Toolbox of MATLAB, such as the FitPercentage, the Akaike’s Information Criterion (AIC), used to avoid overfitting, and the Bayesian Information Criterion (BIC), which penalizes the number of parameters used in the model. In Table 5.2 it is shown a comparison of the fitting that gives the estimated model when using the linear models of Transfer Functions and State Space. In the comparison, only Transfer Functions with relative degree 0 are considered (with the same number of poles and zeros), whose behavior is compared with the equivalent State Space model with the same number of states. p. 52 Report As it can be seen, the more poles/zeros and states used, the better the fitting and the lower the AIC/BIC, so we are interested in looking for big values of poles/zeros and states for the model fitting. Table 5.2: Comparison of the fitting of each estimated model. Parametric model identified % Fit "Next" %Fit "Top" AIC BIC (a) Transfer Function np= nz=8 91.57 % 85.48 % -1.0609e6 -1.0606e6 (b) State space nx=8 states 80.91 % 69.32 % -1.0116e6 -1.0114e6 (c) Transfer Function np=nz=16 99.11 % 96.88 % -1.1662e6 -1.1657e6 (d) State Space nx=16 states 99.03 % 96.18 % -1.1557e6 -1.1554e6 (e) Transfer Function np=nz=19 99.43 % 96.13 % -1.1706e6 -1.1700e6 (f) State Space nx=19 states 99.54 % 97.09 % -1.1929e6 -1.1925e6 Apart from looking at the fitting performance measures, it is also of interest to check the shape of the Bode diagram of the identified models, and compare it with the original one to see how the model differs from the original response. In Figure 5.7 this comparison can be seen, where the same identified models of Table 5.2 are present. In the different subfigures it can be seen that the peaks and valleys that have bigger amplitudes than the others are the ones that are correctly identified even with a small number of parameters, as they have a stronger effect on the frequency response of the system. This can be clearly seen in the first subfigures (a) and (b), where 8 poles and 8 zeros are used to identify the model of ANSYS. The mismatch on the frequency responses is produced due to the poor number of parameters, meaning that the dynamic behavior of the structure cannot be captured just with an 8th order model. Checking the shape of the identified model it can be seen that the central modes with small gain that are in the range of 60-200 Hz are not identified, as it considers a flat frequency response. The approach done is to obtain an identified system as close as possible to the ANSYS model, so it is needed to increase its number of poles and zeros to fit even the ones with a very small gain. Therefore, doubling the number of parameters, the fitting of the identified model is presented in the subfigures (c) and (d), which have almost a perfect fitting, ranging from 99.1% of fit for the "Next" strain gauge to the 99.03% of fit for the "Top". However, the shape of this fitting is not exactly the same as the modeled one, as some peaks of the identified system do not have a correct magnitude. Finally, slightly increasing the number of parameters to 19 it is obtained the fitting presented in (e) and (f). Even though the fitting using a transfer function model is even worse than the fitting using 16 parameters, the fitting using an state space model is almost perfect, having the highest percentage of fit for both strain gauges response (99.54% for "Next" and 97.09 % for "Top"). From this point, increasing the number of parameters only increases the value of the AIC/BIC measure, without having a significant impact to the fitting, so it has been decided to stop increasing the number of parameters and end up choosing 19 parameters. Thus, the final parametric model chosen for representing the system is the one obtained using the State Space representation with 19 states, which corresponds to the fitting of Figure 5.7(f). Modelling and control of the vibration of a rotatory disk p. 53 (a) Transfer Function np= nz=8 (b) State Space nx=8 states (c) Transfer Function np= nz=16 (d) State Space nx=16 states (e) Transfer Function np= nz=19 (f) State Space nx=19 states Figure 5.7: Comparison of the Magnitude (dB) of the identification in air of the "Next" and "Top" strain gauges using different models. p. 54 Report Once the full model with 19 states is obtained, it is needed to perform a model order reduction to keep only the most dominant states and avoid having a big order plant model, that will yield to a big order controller. To do so, several commands in MATLAB exist to reduce the order of a model, such as the balred command, which computes a low-order approximation of the input model. Taking in mind the shape of the identified system of 19 states, it is obtained a reduced order model with 8 states, presenting the frequency response shown in Figure 5.8. As it can be seen, with only 8 states the model is simplified a lot, but it maintains the most relevant poles and zeros, which are the ones with the biggest amplitude change in the Bode diagram. Figure 5.8: Bode diagram of the identified and reduced order model of the disk in air. The MATLAB code to directly obtain the identified model from the tabular data of ANSYS is presented in Script 1. On it, it is shown how the idfrd object is created from the frequency response data imported to MATLAB, where both outputs of the "Next" and "Top" strain gauges are considered in the same Single-Input Multiple-Output (SIMO) system. Then, the final model is estimated using an state space model with 19 states, which is labeled and saved as sys_ss.mat. Finally, it is reduced to an state space model of order 8 using the aforementioned balred command, labeling and saving the system as rsys.mat. Modelling and control of the vibration of a rotatory disk p. 55 Script 1: MATLAB code for model extraction, identification and reduction. 1clc 2clear 3 4% Import the files 5next_100_damp=readtable("SG_next_aire_damping.xlsx",'PreserveVariableNames',true); 6top_100_damp=readtable("SG_top_aire_damping.xlsx",'PreserveVariableNames',true); 7 8% Extract info by variable name 9FREQ_next=next_100_damp.("Frequency [Hz]"); 10 AMP_next=next_100_damp.("Amplitude [m/m]"); 11 PHASE_next=next_100_damp.("Phase Angle [degree]"); 12 13 FREQ_top_100_damp=top_100_damp.("Frequency [Hz]"); 14 AMP_top_100_damp=top_100_damp.("Amplitude [m/m]"); 15 PHASE_top_100_damp=top_100_damp.("Phase Angle [degree]"); 16 17 % Create 2 idfrd models ("sys_next" and "sys_top") 18 19 resp_next=(AMP_next/100).*exp(1i*PHASE_next*pi/180); %Complex value 20 sys_next=idfrd(resp_next,FREQ_next,'FrequencyUnit','Hz'); 21 22 resp_top=(AMP_top_100_damp/100).*exp(1i*PHASE_top_100_damp*pi/180); 23 sys_top=idfrd(resp_top,FREQ_top_100_damp,'FrequencyUnit','Hz'); 24 25 % SIMO with 1 input, 2 outputs ("resp_100") 26 27 Hresp=zeros(2,1,length(FREQ_next)); 28 Hresp(1,1,:)=resp_next; 29 Hresp(2,1,:)=resp_top; 30 resp_100=idfrd(Hresp,FREQ_next,0,'FrequencyUnit','Hz','InputName',... 31 'Voltage','InputUnit','V','OutputName',{'Next','Top'}, ... 32 'OutputUnit',{'m/m','m/m'}); 33 34 % Bode of the SIMO system 35 bode(resp_100,'r'); 36 37 % Identification of the system 38 opt=ssestOptions('EnforceStability',true); 39 sys_ss=ssest(resp_100,19,opt); 40 compare(resp_100,sys_ss) 41 42 % Reduction of the order 43 opt = balredOptions('StateElimMethod','MatchDC'); 44 rsys=balred(sys_ss,8, opt); % Reduce the order selecting 8 states 45 bo = bodeoptions; 46 bo.PhaseMatching = 'on'; 47 bode(rsys,bo) 48 hold on 49 bode(sys_ss,bo) 50 51 % Save the identified and reduced model 52 save('rsys.mat','rsys'); 53 save('sys_ss.mat','sys_ss'); p. 62 Report As it is shown in the table, as long as we increase the number of parameters used, the fit percentage increases and the AIC/BIC decrease, which indicates that the identified model is closer to the frequency response data of the idfrd object. However, there’s a point in which increasing the number of parameters only induces complexity into the identified model and the fitting starts decreasing, as well as increasing the AIC/BIC. As it can be clearly seen in the tabular data, this inflection point occurs in the identified model with 19 parameters, number from which increasing the parameters only decrease the fit percentage of all the outputs and slightly increase the AIC/BIC values. Hence, as the interest relies on the production of an identified model as complete and reliable as possible, the 19 state space model is used. The fitting of each parametric model presented on the table is shown in Figure 5.14, where the fitting of both "Next" and "Top" strain gauges are compared with the original system. (a) Identified model with nx=8 parameters. (b) Identified model with nx=16 parameters. (c) Identified model with nx=19 parameters. (d) Identified model with nx=20 parameters. Figure 5.14: Comparison of the Magnitude (dB) of the identification in water of the "Next" and "Top" strain gauges using different models. We can see that with only 8 parameters (subfigure a) the model is able to represent just some of the peaks and valleys produced at the natural frequencies, but is not enough for its correct identification. Doubling the number of parameters to 16 (subfigure b), the fitting improves, Modelling and control of the vibration of a rotatory disk p. 63 but some zones of the frequency response are still unable to correctly follow the same gain and shape of the original model. Finally, with 19 parameters (subfigure c) the model is perfectly fit, except for one tiny natural frequency between 60 Hz and 70 Hz. Trying to increase the number of parameters to 20 (subfigure d) may be able to fit that zone but, in exchange, the accuracy of the overall system decreases. Therefore, the final number of parameters used to identify the plant model in water is 19 parameters. The fit percentage of the state space model with 19 states is bigger than the transfer function model with 19 poles and 19 zeros, so it is finally identified the plant model using the state space model presented in Figure 5.14 (c). The obtained plant with 19 states represents the behavior of the modeled system, but it is not suitable for designing a controller, as it will yield to a controller of order 19, which is too big. Therefore, once obtained the reliable model of 19 states, a model order reduction is performed to keep only the most dominant behavior of the system and get rid of the small gain peaks of the frequency response that won’t be significant in the vibration of the disk inside the water. Therefore, using the same order reduction method as presented in Script 1, the reduced order model is obtained. In this case, more natural frequencies appear in the range of 10-326 Hz than in the air, so it is needed more states in the reduced order model to capture the dominant behavior at that resonant frequencies. Therefore, 12 states are used to describe the reduced order model of the system in water, which is compared with the identified system of 19 states in Figure 5.15. Figure 5.15: Bode diagram of the identified and reduced order model of the disk in water. p. 64 Report Finally, the plant used for the design of the controller of the disk under water is defined by the state space representation of equation (5.13), where the state space matrices are the following ones. Notice that in order to improve the representation of the 12-dimensional matrices, the values of the matrices have been rounded to the nearest integer. However, this is only a representation rounding, as the matrices used in MATLAB include all the corresponding decimal numbers. ˙x=A2x+B2u y=C2x+D2u(5.13) A2=                      −57 1345 −227 −2185 35 275 31 −42 −906 −18 −53 124 −1434 33 −1320 33 −585 −201 −89 163 52 −80 44 −105 98 1234 −58 −2075 318 127 −102 −209 −879 130 19 115 2389 15 2226 −23 791 180 269 −92 −60 42 −108 −51 −21 553 −339 −758 −9−85 72 193 −428 105 −35 74 −460 162 −66 −174 22 95 −381 −337 −159 −79 112 −37 −142 22 85 −184 −49 448 10 −164 −770 −441 −163 193 14 −121 180 157 67 −36 −27 −230 904 −196 83 −214 759 −78 896 −17 457 240 695 −1161 16 42 −148 83 182 77 −160 6 −23 −180 210 572 −69 169 −172 382 53 −12 11 −30 −145 −4 201 103 180 −95 −4−1331 −73 107 −22 281 −30 46 −164 −19 −41 −466 899 −192                      B2=                      21.55 70.59 12.97 −117.18 15.03 20.57 4.58 −10.73 −65.38 −3.13 −8.19 73.42                      ×10−4 C2="−615 81 −559 −290 −1158 −468 −219 −94 −187 −129 181 29 −600 137 −626 −83 −36 27 −28 160 −87 −38 108 9 #×10−5 D2="−4.4661 −2.4816 #×10−9 Modelling and control of the vibration of a rotatory disk p. 65 5.3 Comparison of the plant of the disk in air and water First, it is important to point out the differences between the frequency responses obtained from both simulations in air and water, which are the models for the identified plants. In Table 5.5 it is summarized the different natural frequencies obtained from the harmonic analysis of the system in air and water for each strain gauge. As commented before, it can be seen that the natural frequencies when the disk is submerged in water appear in lower frequencies than in air, which is due to the added mass effect induced by the tank of water. Table 5.5: Comparison of the natural frequencies obtained from the "Next" and "Top" strain gauges in air and in water. Type of mode shape Natural frequency "Next" Natural frequency "Top" Harmonic air (Hz) Harmonic water (Hz) Harmonic air (Hz) Harmonic water (Hz) 0NC 49.952 14.982 49.952 14.982 1ND 43.167 16.382 43.167 16.382 2ND 64.797 27.512 64.797 27.512 3ND 133.27 62.5 133.27 63.5 4ND 232.9 118.0 233.1 118.2 1NC 305.8 126.772 305.8 126.772 1NC-1ND 330.606 135.118 330.606 135.118 In order to quantify the added mass effect induced by the water, it is computed an adimensional value called Frequency Reduction Ratio (FRR), which is obtained by comparing the natural frequencies in air and water, according to equation (5.14). FRR (%) =fair −fwater fair ×100 (5.14) The computed FRR for the first mode shapes of the disk take values between 49.29 % and 70.01 %, as presented in Table 5.6. It can be seen that the FRR (%) decreases as the number of nodal diameters (ND) increase, which is produced because the deformation of the mode shapes move less water than before. This produces that the value of the added mass effect will be lower, having a smaller affectation to the frequency. The same effect happens when the number of nodal circles (NC) increase from 0NC to 1NC, reducing the vertical displacement of the disk and, hence, decreasing the quantity of water displaced. Table 5.6: Computation of the Frequency Reduction Ratio (FRR) for the first mode shapes of the structure. Type of mode shape FRR "Next" (%) FRR "Top" (%) 0NC 70.01 % 70.01 % 1ND 62.05 % 62.05 % 2ND 57.54 % 57.54 % 3ND 53.10 % 52.35 % 4ND 49.33 % 49.29 % 1NC 58.54 % 58.54 % 1NC-1ND 59.13 % 59.13 % p. 66 Report Additionally, the frequency response of the system in air and water can be compared in the same bode diagram, as presented in Figure 5.16. On it, it can be seen the effect of the reduction of the natural frequencies, which produces that the frequency response of the plant model in water is concentrated in low frequencies, having more peaks in the range of 10-326 Hz. Moreover, it can be seen that both plants have the same DC gain, which means that the step response of both systems will have the same steady-state value, corresponding to the value of the transfer function evaluated at s= 0, as stated by the Final Value Theorem. Figure 5.16: Comparison of the bode diagram of the plant model of the system in air and in water. Modelling and control of the vibration of a rotatory disk p. 67 6 Experimental validation of the model In order to validate experimentally the identified plant models obtained from the ANSYS simulation, we performed different experiments on the test rig mounted on the laboratory to extract the natural frequencies of the system and compare them with the ones obtained from the bode diagram of the identified models. To do it, different sensors and piezoelectric patches are placed on different positions of the disk. Among all the sensors and piezoelectric patches placed, we can see the "Top" and "Next" strain gauges used in the simulations in Figure 6.1a and Figure 6.1b, together with the piezoelectric patch used. Additionally to the "Top" and "Next" strain gauges, two more gauges are placed in the top and bottom of the disk, as it can be seen in the image. These two extra gauges are placed near to the shaft of the disk, separated 45◦from the "Next" and "Top" gauges. The two additional gauges are are also monitorized in the laboratory, and are labeled as "Next_2" and "Top_2". (a) Top view of the disk with the sensors. (b) Bottom view of the disk with the sensors. Figure 6.1: Zoom in of the sensors placed on the disk. An important thing to mention is that the piezoelectric patch used to validate the numerical model is a different one from the patch simulated in the previous section, corresponding to the "P-876.SP1" instead of the "P-876.A12". This change was due to the disponibility of the patches on the laboratory, and it is outside of the scope of the project, as it depends on the laboratory physical setup. However, both piezoelectric patches have the same operation range of +400 V to -100 V and they are made of the same PIC255 material, so they only differ on the size (the SP1 mounted on the disk is smaller than the A12). Therefore, it has been changed the solid representing the piezoelectric patch in the simulation to adapt to the new geometry, obtaining new results different from the ones presented on the previous section. p. 68 Report Additionally, two accelerometers are used to measure the vibration of the structure and obtain a more reliable result of the position of the natural frequencies of the disk. As shown in Figure 6.2, they are placed on the tip of the disk, where the displacement is bigger. Figure 6.2: Top view of the disk, with the accelerometers marked. The sensors and the piezoelectric patch have to be connected to the corresponding National Instruments (NI) modules of the CompactDAQ "cDAQ9185" chassis, used to receive and transfer the measured data to the computer, presented in Figure 6.3a. Concretely, 4 modules are connected to it, being the first module used to obtain the acceleration from the accelerometers. The second module is an input-voltage module, used to monitor the voltage sent to the patch, whereas the third one is an output-voltage module, used to generate a voltage that ranges from 8 V to -2 V. This voltage is then amplified by a factor of 50 using a piezo amplifier, which also serves as a driver to control the voltage sent and avoid burning the patch with an excessive voltage out of its operational limits. Finally, the fourth module is used to measure the strain from the strain gauges placed on the disk ("Next" and "Top", among others). Once connected the modules to the cDAQ, we have to program its acquisition by using an Ethernet cable and a program called NI MAX (Measurement and Automation Explorer). Using that program, we can define the different characteristics of the input sensors, such as the measurement range, the sensitivity and electrical connections. In Figure 6.3b it can be seen the configuration used for the strain gauges, whose strain configuration is set to a quarter bridge because the Wheatstone bridge that is used to measure the strain is partially implemented inside the NI module [25]. Additionally, the gauge factor provided by the manufacturers is introduced into the program, which is internally used to compensate the Wheatstone bridge of the measurement. Regarding the lead resistance, it represents the resistance induced from the wires. Hence, from the datasheet of the strain gauges, we obtain that the configuration 7/0.12 of the gauges has a total resistance of 0.44 Ω/m, with a cross section diameter of the wires of 0.08 m2. As the gauges have 3 meters of wire, the total lead resistance introduced is 1.32 Ω. Modelling and control of the vibration of a rotatory disk p. 69 (a) CompactDAQ and modules used to connect the sensors and send voltage to the patch. (b) NI MAX programming of the different sensors used. Figure 6.3: Hardware and software needed to obtain experimentally measures from the different sensors. Finally, the strain gauges have to be calibrated to compensate the small initial deformation that the gauge should have when installed, as it is impossible to install it perfectly parallel to the disk surface. This is done by sending an initial voltage to the Wheatstone bridge, forcing a strain measured of 0 when the disk is completely still. Once programmed the sensors, the numerical validation is done by sending a chirp signal of maximum amplitude to the piezoelectric patch, which is a signal in which the frequency increases linearly with time. This chirp is generated programming the signal using LabVIEW, and sent to the patch by the corresponding output-voltage module connected to the compactDAQ. By exciting the piezoelectric patch with this frequency-varying signal, it can be determined at which frequency the structure starts vibrating because it will correspond to the frequency at which it will increase the strain and acceleration measured, producing different peaks on the temporal signal measured. In Figure 6.4 it is presented a slice of 2 s of the chirp signal monitorized by the LabVIEW program, as well as the acceleration measured exciting the patch with the full chirp of 120 s, where each time that increases its amplitude its due to the vibration of the structure produced by a natural frequency. (a) Slice of 2 s of the chirp signal monitorized. (b) Acceleration signal measured in 120 s. Figure 6.4: Example of the chirp and acceleration measured by the sensors using LabVIEW. p. 70 Report 6.1 Validation of the disk in air The chirp signal has been tested when the structure is in air and when it is submerged at different heights inside the tank of water, in order to validate both models. For both cases, the acceleration and strain have been measured using an acquisition program of LabVIEW that stores the data in a Technical Data Management Streaming (TDMS) file. Then, in order to identify the natural frequencies, another LabVIEW program has been implemented to read the previous TDMS file, with all the collected data, and perform the FFT spectrum of the signals obtained to check their peaks of amplitude. The acquisition rate of the measured signals is set to 5000 Hz, which is the same as the rate of the generated chirp signal. For defining the duration of the chirp, which will determine the ramping speed of the frequency, it has been tested different chirps of 15, 30, 60 and 120 s. We saw that the longer the chirp, the better the FFT results were obtained, so the approach done experimentally was to perform 5 different chirps of 120 s to cover the range of 0-550 Hz, obtaining data in the range of 0-150 Hz, 100-250 Hz, 200-350 Hz, 300-450 Hz and 400-550 Hz. Therefore, a first signal ranging from 0-150 Hz in 120 s was performed, obtaining the FFT shown in Figure 6.5 for the case of the structure in air. Figure 6.5: FFT of the acceleration and strain obtained from the chirp 0-150 Hz in air. As the strain and the acceleration have different orders of magnitude, two different axis are used to plot both measures in the same plot, corresponding to the strain (m/m) in the left axis and the acceleration (m/s2) in the right axis. In a first sight it can be clearly seen that whereas the acceleration signal (red and green lines) detect the changes of amplitude at different frequencies, the strain gauges (purple and black lines) detect a peak value at 50 Hz and the rest of the spectrum is very noisy. This happens because strain gauges are like antennas, meaning that they are very sensitive to the electrical noise of 50 Hz. However, if we zoom in around a natural frequency, as shown in Figure 6.6, we can see that the strain clearly increases when a particular mode shape is excited. The peaks of strain occur in the same frequency as the peaks of the acceleration, so we can ensure that a natural frequency appears at that point. In the frequency range of the figure, three different natural frequencies can be seen, corresponding to the three peaks of the computed FFT. The first two peaks occurring at 35.48 Hz and 37.22 Hz correspond to the double modes of the 1ND mode shape, as both peaks of strain appear in different frequencies and are rotated a certain angle. This rotation of the Modelling and control of the vibration of a rotatory disk p. 71 nodal diameter produces that depending on the position of the sensors and the mode shape, one sensor may be able to detect better one frequency or the double mode of it. The third mode shape that appears in 38.94 Hz produces that the two strain gauges measure the same peak of strain, meaning that it is a circumferential mode that induces the same strain in circles. Concretely, it is identified the 0NC mode shape in that frequency of 38.94 Hz. Following the same approach as in the modal analysis, these double modes are going to be simplified by taking the average value of both frequencies and considering them as a unique mode occurring at that frequency. Hence, the previous double modes are considered to be appearing at an average frequency of 36.35 Hz, corresponding to a unique 1ND mode shape. Figure 6.6: Zoom in of the first mode shapes obtained from the chirp 0-150 Hz in air. Comparing the FFT obtained from the accelerometers with the FFT obtained using the strain gauges, it can be seen that with the current level of excitation the peaks of strain obtained are very small compared with the acceleration peaks. This fact implies that the detection of the mode shapes using these strain gauges is not going to be as accurate as using the accelerometers to determine that experimental modal analysis. One way to increase the peaks obtained by the gauges could be to use a more powerful piezoelectric patch, to induce more strain to the disk and increase the amplitude of the mode shape. On the other hand, to detect better these small changes in strain, we could use different strain gauges with more resolution. Due to that limitations of the strain gauges, the extraction of the natural frequencies to validate the numerical model is going to be used based on the acceleration signal, as with the current level of excitation they give more accurate results than the gauges. Moreover, they are less affected by the electrical noise and they are able to detect all the mode shapes of the structure even though they do not produce strain in the direction of the gauges. Computing the FFT of the results obtained using the other chirps in the frequency range of 0- 550 Hz, and checking the frequency where a maximum occurs, we can obtain the experimental natural frequencies of the disk in air, neglecting the peaks of electrical noise at 50 Hz and multiples. Finally, the collected experimental values are compared with the ones obtained in the numerical simulation with the new geometry, as presented in Table 6.1. As stated before, these simulated values differ a little bit from the ones presented in the previous section, but can be used anyway to demonstrate the validity of the simulation against the experimental data, as only the geometry has been changed in the simulation. p. 78 Report over, it is categorized as optimal because it tries to reach the performance specifications while reducing the amount of work done by the control system, process which is done by optimizing a cost function J(u). This type of controller is implemented using a full state feedback gain, so it is needed to obtain all the states of the system. However, only the output of the "Next" strain gauge is considered to be obtained from the system, so it is needed to design an observer together with the LQR in order to estimate all the states of the system and implement it in the laboratory, obtaining an LQG controller. On the other hand, the second control technique implemented is a more manual control design, presented in Subsection 7.3, which consists on directly cancelling the first poles and zeros of the frequency response of the system. This is done by adding in the direct control chain the zeros and poles corresponding to the first poles and zeros of the plant, respectively. Together with this control system, it has been added a first order filter to attenuate the effect of the high-frequency vibrations that are not directly cancelled. 7.1 Controllability and observability of the system The reduced order state space model for the simulation of the disk in air has 8 states, 2 outputs, and 1 input, so the Amatrix is an 8×8state matrix, Bis an 8×1input matrix, Cis a 2×8 output matrix and Dis a 2×1feedforward matrix. To check if the system is controllable, the controllability matrix Ris computed using the following equation (7.1), where nis the number of states, and the matrices Aand Bare the corresponding state space matrices of the model (5.10). R=hB AB A2B A3B . . . An−1Bi(7.1) Once the 8×8controllability matrix Ris computed, the controllability of the system is checked by looking at its rank, which needs to be maximum, meaning that the rank(R) = n= 8. On the other hand, to check if the system is observable, the observability matrix Ois computed using equation (7.2), where Aand Care the matrices of the same state space model. O=           C CA CA2 CA3 . . . CAn−1           (7.2) In this case the system is observable if the observability matrix has full rank, which implies that the rank(O) = n= 8. The MATLAB code to load the identified state space models and compute the ranks of those controllability and observability matrices is presented in Script 2. With it, it can be checked the controllability and observability of the reduced order model, which returns that effectively the reduced order system modeled from the ANSYS simulation in air is controllable and observable. Modelling and control of the vibration of a rotatory disk p. 79 Script 2: MATLAB code to check controllability and observability matrices. 1% Load the saved reduced order system and full−order system: 2load('rsys.mat'); 3load('sys_ss.mat'); 4 5% Obtain state space matrices from reduced order model: 6A=rsys.A; 7B=rsys.B; 8C=rsys.C; 9D=rsys.D; 10 11 R=[B,A*B,A^2*B,A^3*B,A^4*B,A^5*B,A^6*B,A^7*B]; % Controllability matrix 12 disp('Rank R: ') 13 disp(rank(R,10^−6)) % tolerance of 10^−6 14 o=[C;C*A;C*A^2;C*A^3;C*A^4;C*A^5;C*A^6;C*A^7]; % Observability matrix 15 disp('Rank O: ') 16 disp(rank(o,10^−6)) % tolerance of 10^−6 In the case of the system submerged in water, the computation of both controllability and observability matrices is done considering the new rank of the reduced order model in water, which is n= 12. Then, the same script can be used to import their corresponding plant models and check that, effectively, the rank of those matrices is also maximum. All in all, it has been demonstrated that both plant models of the system in air and water are controllable and observable, so the different control techniques can be applied to those systems to reduce the vibrations and the strain of the disk. 7.2 Control approach 1: LQR controller + observer The Linear Quadratic Regulator (LQR) is a type of optimal control that is based on state space representation, and it is used to obtain a full state feedback gain by choosing closed-loop characteristics that are important for the system. Specifically, it takes into account how well the system performs and how much effort does it take to get that performance. This is done by defining a cost function J(u)that adds up the weighted sum of performance and effort over all time and tries to minimize it. J(u) = Z∞ 0xTQx +uTRu +  2xTNudt (7.3) As presented in equation (7.3), the states xand the control action uare multiplied by two matrices Qand R, being Qa symmetric positive semidefinite matrix (such that Q≥0) and R a symmetric positive definite matrix (R > 0). Both Qand Rmatrices represent the weights used to penalize the states xand the input of the system u, giving more priority either to the minimization of the state tracking error or to the minimization of the control action needed, respectively. Additionally, an Nmatrix can be included to penalize cross products of the input uand the states x, but in order not to add extra parameters to the optimization problem, it is just going to be set Nto zero and focus only on Qand R. p. 80 Report In our case, the desired output of the system yto be controlled is the strain read by the "Next" strain gauge, which has an order of magnitude of around 5×10−6m/m. As a consequence, the states of the system xare also of that small order. On the contrary, the input of the system uis the voltage of the piezoelectric patch, which ranges from a value of +400 V to -100 V. Therefore, the weights Qand Rare first used to normalize that two magnitudes to be able to correctly weight each other, ending up with the weights presented in equation (7.4), where the same priority is given to the minimization of the states and the control action. Q=1 5×10−6−0= 2 ×105 R=1 400 −(−100) = 2 ×10−3(7.4) In order to give more priority to the minimization of the strain, it has to be taken into account that when the Rmatrix is kept constant and the value of Qincreases, the system will evolve faster to the reference, obtaining bigger control actions uand having an aggressive behavior. On the contrary, if the weight of the control action Ris increased, keeping Qconstant, the control action uwill take smaller values in order to save energy, but the system will take more time to evolve towards the reference, having a conservative behavior. After searching for a trade-off between both performances, the weight for the tracking of the states has been set to be ten times bigger than the weight of the control action, ending up with the final weights of Q= 2 ×106 and R= 2 ×10−3. This tuning is shown graphically after simulating the results of the system using that weights, as presented in the further Figure 7.5. With the LQR design, the optimal control action obtained corresponds to the state-feedback control law u=−Kx, which depends on an optimal gain matrix K, obtained as a result of the minimization of the quadratic cost function J(u), subject to the dynamics of the system: ˙x=Ax +Bu y=Cx +Du (7.5) The way that MATLAB solves the LQR problem is by deriving Kfrom the Algebraic Riccati’s equation solution Spresented in equation (7.6), which must satisfy that S≥0, and then use equation (7.7) to obtain the gain K. Notice that the term Nhas been introduced in the equations but is not considered in this case, as commented before. ATS+SA −(SB +  N)R−1(BTS+  NT) + Q= 0 (7.6) K=R−1(BTS+  NT)(7.7) Finally, the LQR is implemented closing the loop and checking for the stability of the closed-loop system, looking at the poles of Acl = (A−BK), values which are also returned by the MATLAB command lqr(A,B,Q,R). By default, this controller makes the states of the closed-loop system tend to 0, as the dynamics of the system are defined by ˙x= (A−BK)x. In order to define a non-zero output for the controlled system, it has to be added an additional input reference rto track, as shown in the control scheme of Figure 7.1. This input reference rrepresents the desired "safe" strain that the structure has to attain to stop vibrating, so it has to be tracked to accomplish the control objective. Modelling and control of the vibration of a rotatory disk p. 81 With this extra input, the control law ends up being u=−Kx +v, where v=k∗r. This precompensation gain k∗is added after the input reference rand before the closed-loop of the system, and it is computed as the inverse of the gain of the closed loop system, making the output yto converge towards rinstead of towards 0, as the dynamics of the system are now defined by ˙x= (A−BK)x+Bv. Figure 7.1: Scheme of the LQR controller implemented in Simulink. 7.2.1 LQR controller on the plant model of the disk in air The effect that closing the loop has on the frequency response of the system in air it is shown in the Figure 7.2, where the Magnitude and Phase of the frequency response of the "Next" strain gauge are compared with the original plant model without the controller. By visual inspection, it can be seen that this LQR controller reduces the magnitude of all the peaks of strain at the natural frequencies, attenuating them, while it keeps the original magnitude of the rest of the bode diagram. This produces that the structure is not going to vibrate as much as before at the natural frequencies, as the amplitude of the vibration is reduced. Figure 7.2: Effect of the LQR controller in the frequency response of "Next" in air. p. 82 Report The effect of this attenuation of the frequency response can be seen in the temporal domain by computing the step response of the closed-loop system, which corresponds to the temporal evolution of the strain read by "Next" strain gauge when the input voltage of the piezoelectric patch changes from 0 V to 1 V. This step response is presented in Figure 7.3, in which it can be seen that the closed-loop system with the LQR controller converges faster to the steady state value of -4×10−8m/m. Moreover, the high-frequency oscillations that were present in the system without the controller are also eliminated with the designed controller. Figure 7.3: Effect of the LQR controller to the step response of the the "Next" strain gauge in air. The tracking of the "safe" strain reference value is checked using the previous Simulink model of Figure 7.1, in which the reference value rcorresponds to a step of 1×10−6m/m at time 20 s. Introducing the state feedback gain obtained from the LQR design, it can be simulated the evolution of the strain measured by the "Next" and the "Top" strain gauges towards that "safe" strain value, which is presented in Figure 7.4. Figure 7.4: Tracking of a "safe" strain of 1×10−6using LQR without observer in air. Modelling and control of the vibration of a rotatory disk p. 83 In the figure it can be seen how the strain measured by the "Next" strain gauge converges towards the reference value in approximately 0.2 s, while the "Top" strain gauge remains stable around the corresponding strain of that part of the disk. The control action needed to achieve that behavior is presented in the right hand side of the figure, which is indeed between the operational limits of the piezoelectric patch (+400 V to -100 V), being a feasible value. The influence of the tuned weights Qand Ron the temporal evolution of the outputs and the control actions is presented in Figure 7.5. On it, it can be clearly seen that, once the matrices are normalized using the values of equation (7.4), they can be used to prioritize the minimization of the tracking error (increasing the weight of Q) or the minimization of the control action (increasing the weight of R). In this case, the final weights used correspond to the blue line, which is to set the minimization of the tracking error to be 10 times more important than the minimization of the control action, as is more important to minimize the oscillations of strain. Figure 7.5: Tuning of the Q and R weights using the temporal response of the closed-loop system in air. p. 84 Report 7.2.2 LQR controller on the plant model of the disk in water On the other hand, the effect that the LQR controller has on the frequency response of the plant model of the system submerged in water is presented in Figure 7.6, where the Magnitude and Phase of the frequency response of the "Next" strain gauge are also compared with the original plant model without the controller. In this case, it is also clear that this controller attenuates the vibration of the disk under the water by reducing the peaks of strain at the natural frequencies of the disk. Figure 7.6: Effect of the LQR controller to the frequency response of the "Next" strain gauge output of the system under the water. This attenuation of the frequency response can also be seen in the temporal domain by computing the step response of the closed-loop system, shown in Figure 7.7. In the figure, it can be seen how the controlled system converges faster towards the steady-state value without all the high frequency values that were present in the step response of the system without the controller. Figure 7.7: Effect of the LQR controller on the step response of the disk under water. Modelling and control of the vibration of a rotatory disk p. 85 In the case of simulating the system submerged under the water, the same Simulink model presented in Figure 7.1 is used but changing the plant and the optimal gain Kobtained, ending up with the tracking of the reference 1×10−6m/m shown in Figure 7.8. With this controller and with the submerged system in water, the "Next" strain gauge takes longer to converge towards the reference, taking approximately 1.2 s. However, a bigger overshoot value is present on the response, which rises the strain up to 1.8×10−6m/m before converging to the safe value set of 1×10−6m/m. Regarding the voltage needed, it is inside the operating range of the piezoelectric patch. Figure 7.8: Tracking of a "safe" strain of 1×10−6m/m using LQR without observer in water. 7.2.3 Implementation of the Kalman filter observer Even though in a simulated model it can be implemented the full-state feedback of the designed LQR controller, it is not possible to implement it in the laboratory because all the states of the system cannot be directly measured with the sensors available. However, this problem can be solved using an observer to estimate those state values from the measured output of the disk, which is the "Next" strain gauge. In the previous Subsection 7.1 it has been demonstrated that both plant models are observable, so it will be designed an observer together with the previous LQR controller to estimate the states of the system ˆxusing the outputs of the plant y. This Luenberger observer is parameterized using the matrix Lof static gains, which applies a correction to the estimated states ˆ ˙x according to the output measured y, as it is shown in equation (7.8). ˆ ˙x=Aˆx+Bu +L(y−Cˆx)(7.8) p. 86 Report The design of an observer is dual to the state feedback design, but instead of considering the A and Bmatrices to obtain the feedback gain K, it is considered A0and C0matrices to obtain the observer gain L. In this case, the observer is designed using a Kalman filter, which considers that the state space model of the system is augmented with two additional inputs, as presented in equation (7.9). These two inputs model the disturbance entering in the control action w, which has a direct affectation to the system, and the noise vof the sensors, which is a measurement noise that affects the output. ˙x=Ax +Bu +Bww y=Cx +Du +v(7.9) These system disturbance w(t)and measurement noise v(t)are defined as uncorrelated Gaussian noise processes with zero-mean, where the associated disturbance and measurement noise covariance matrices Wand Vare defined as presented in equation (7.10). E{wwT}=W≥0, E{vvT}=V≥0, E{wvT}= 0.(7.10) The Kalman filter will then obtain an optimal estimation of the state ˆx, which minimizes the expected value E{(x−ˆx)T(x−ˆx)}. Using these estimated states instead of the real unknown states, the control action is computed with the full-state feedback gain obtained with the LQR approach. The coupling between the observer (state estimator) and a state feedback gain K (LQR) is programmed using a regulator, as it is shown in the scheme of Figure 7.9. Figure 7.9: Scheme of the regulator with an observer and a feedback gain. Image from: [3] The coupling between the Kalman filter and the LQR optimal control is known as the Linear Quadratic Gaussian (LQG) optimal control. In MATLAB, it is computed using the lqgreg command, which allows to merge the optimal state feedback gain Kwith the parameterized observer gain L. Then, the resulting regulator Klqg is directly used as an output feedback controller, connected to the measured output yof the system, as shown in the scheme of the regulator. Modelling and control of the vibration of a rotatory disk p. 87 7.2.4 LQG controller (LQR controller + Kalman filter observer) on the plant model of the disk in air The implementation of this equivalent regulator in MATLAB ends up modifying the frequency response obtained with the previous LQR control, as the Kalman filter changes the dynamics of the closed-loop system. For the system in air, the effect of this combined LQG control is shown in Figure 7.10, where the bode plot of the frequency response of the LQR controller and the original system are also represented. As it is shown, the new LQG controller (LQR + observer) changes a little bit the bode diagram of the closed-loop system, increasing the gain of the system in some zones. However, the important issue is that the peaks of the natural frequencies are still attenuated, reducing the vibration of the structure. Figure 7.10: Comparison between effect LQR and LQG on the frequency response in air. In order to simulate the performance of this controller, a new Simulink program has been used with different noise and disturbances, which is presented in Figure 7.11. The noise vhas been added to the measurement y1of the "Next" strain gauge, representing the measurement errors of the sensor, and the disturbance whas been added entering in the control action, affecting the plant dynamics. It can be seen how the LQG regulator Klqg is directly considered as an output feedback controller, which takes the signal y1read by the "Next" strain gauge to estimate all the states ˆxof the plant, and uses the full-state feedback gain to compute a new control action u. p. 94 Report Figure 7.19: Effect of the cancellation of poles and zeros to the frequency response of the the "Next" strain gauge output under water. The poles and zeros of the reduced model submerged in water are presented in (7.12), where the 1ND mode shape corresponds in this case to the second pole p2, occurring at frequency f2= 16.3811 Hz, whereas the 0NC mode shape corresponds to the first pole p1, being at frequency f1= 14.9811 Hz. p1=−0.9433 ±94.1242i↔f1= 14.9811Hz p2=−1.0297 ±102.92i↔f2= 16.3811Hz z1=−0.9521 ±80.0817i↔f3= 12.7463Hz z2=−1.0030 ±98.4269i↔f4= 15.6659Hz (7.12) Additionally to these controllers, a first order filter is introduced in the control system to attenuate the high frequency components that are unaffected by the cancellation of the first poles and zeros. This is done to avoid increasing the order of the controller with several poles and zeros of the plant, so the approach is not to directly cancel them, but design a low-pass filter to avoid the strain amplification in those high uncancelled frequencies. The filter designed has the structure presented in equation (7.13), where the gain Khas been set to 1 to avoid modifying the closed-loop gain (K= 1) and the wbparameter has been tuned for each scenario. Hfilter(s) = K s wb+ 1 (7.13) Modelling and control of the vibration of a rotatory disk p. 95 The tuning of the wbparameter is done considering that it defines the bandwidth of the filter, meaning that the filter will have a magnitude gain of -3 dB in the frequency wb, and will then attenuate the higher frequencies with a slope of -20 dB/decade. As we want the first zone of the 1ND and 0NC mode shapes to have a constant slope, it is designed considering that the cut-off frequency has to be before the high-frequency components of the new closed-loop system. In the case of air, the first uncancelled natural frequency is at around 300 Hz, with a gain of approximately 10 dB above the DC gain. Therefore, to attenuate this gain of the high-frequency components, we select a cut-off frequency of wb= 1000 rad/s = 159.2Hz. By doing so, the frequency response of the closed-loop system is modified as presented in Figure 7.20, where the gain of the high-frequency components is now below the DC gain. The addition of this filter to the already cancelled system without the first poles and zeros do not affect the stability of the closed loop, it just improves the performance of the closed loop system. Figure 7.20: Comparison of the frequency response of the the "Next" strain gauge in air with the filter to attenuate high frequency components. On the other hand, when the system is submerged in water, the first uncancelled natural frequency appears at around 130 Hz, with a gain of around 10 dB above the DC gain. Therefore, the tuning of the cut-off frequency in this case ends up being wb= 500 rad/s = 79.58 Hz. By doing so, the frequency response of the closed-loop system considering the controller with the filter under water is modified to the one presented in Figure 7.21. In a similar way, this filter attenuates the high frequency components to reduce their gain below the DC gain of the system. p. 96 Report Figure 7.21: Comparison of the frequency response of the the "Next" strain gauge under water with the filter to attenuate high frequency components. The simulation of this controller is done using the Simulink program presented in Figure 7.22, where the reference signal rand the precompensation gain k∗are used to track the "safe" strain signal. Then, the filter Hfilter and the controller itself are added in the direct chain to cancel the corresponding poles and zeros, computing the control action uto input to the piezoelectric patch of the disk. Figure 7.22: Scheme of the Simulink model used to implement the cancellation of poles and zeros controller and the low-pass filter. In this case, the tracking of the "Next" strain gauge is shown in Figure 7.23, where the convergence of the "Next" output towards the "safe" strain reference is very fast. However, the other output, whose first harmonics are not being cancelled, has a convergence towards its corresponding steady state value slower, but it remains stable. Modelling and control of the vibration of a rotatory disk p. 97 Figure 7.23: Tracking of the reference with the cancellation of poles/zeros and the filter in air. As the poles of the system are intrinsic, the controller designed in this second approach will cancel the poles of both "Next" and "Top" strain gauges outputs using the zeros of the controller. However, the zeros of the system depend on the output, so each gauge has different zeros. This produces that the controller will only cancel the zeros of the "Next" strain gauge, which are not the same as the zeros of the "Top" strain gauge. As a consequence, the poles of the controller introduced to cancel the zeros of the "Next" strain gauge will appear in the frequency response of the "Top" strain gauge, as presented in Figure 7.24. Figure 7.24: Effect of the cancellation of poles and zeros of the "Next" strain gauge to the "Top" frequency response in air. Therefore, those poles will produce an extra oscillation of the "Top" strain gauge, whose first peak of magnitude appears in 36.38 Hz, corresponding to the frequency of the cancelled zero of the "Next" strain gauge presented in (7.11). This produces that, by using this controller, the "Top" strain gauge output will present oscillations at that frequency, as presented in the previous p. 98 Report Figure 7.23. Measuring the frequency of that oscillations, we obtain an approximate value of 35.84 Hz, which is very similar to the first peak of the frequency response of 36.38 Hz. When the disk is submerged in water, the convergence of the "Next" strain gauge is a little bit slower, as it can be seen in Figure 7.28. Figure 7.25: Tracking of the reference with the cancellation of poles/zeros and the filter under water. In this case, the "Top" strain gauge oscillates at a lower frequency than in air, ending up converging to the steady state value after approximately 4 s. In a similar way, that oscillations are due to the poles introduced in the controller to cancel the zeros of the "Next" strain gauge, which modify the frequency response of the "Top" strain gauge, as presented in Figure 7.26. Figure 7.26: Effect of the cancellation of poles and zeros of the "Next" strain gauge to the "Top" frequency response under water. Modelling and control of the vibration of a rotatory disk p. 99 The oscillations of the "Top" output of the previous Figure 7.25 have an approximated measured frequency of 12.53 Hz. As expected, those oscillations are due to the first cancelled zero of the "Next" strain gauge at 12.75 Hz, as presented in the previous equation (7.12). 7.3.1 Test the control approach 2 against external disturbances In order to check how these controllers reacts to an external disturbance in air or water, it has been added an impulse disturbance of 10 V in the control action at time t=1 s, as shown in Figure 7.27 for the case of the disk in air. This disturbance produces that the "Next" strain gauge has to converge another time to the reference value, being able to successfully reduce the vibrations of the disk produced by that disturbance. Figure 7.27: Rejection of an impulse disturbance using the cancellation poles/zeros in air. In the case of the simulation under the water, presented in Figure 7.28, the disturbance has been added at t=2 s to leave more time to the "Top" strain gauge to converge, as it oscillates at lower frequencies. In a similar manner, the increase of strain produced by the disturbance of 10 V has also been corrected and the outputs end up converging another time to the "safe" strain value of 1×10−6m/m. When the disturbance is applied in both air and water scenarios, the "Next" strain gauge oscillates more than using the first LQG control approach. This happens because this cancellation of poles and zeros is not designed to minimize the impulse frequency response, as it was the case of the LQG approach of the first controller. However, it is still able to track another time the "safe" strain value once some external factor increases that strain value, so it still works as long as the disturbance is an impulse. p. 100 Report Figure 7.28: Rejection of an impulse disturbance using the cancellation poles/zeros in water. 7.4 Comparison between approaches The main difference between both controllers is their complexity, being the first LQR controller of order 8 or 12, depending on the reduced order model used as the plant of the system. On the other hand, the second controller has order 4 for both air and water scenarios, as it only contains the 2 poles and the 2 zeros needed to cancel the first vibrational modes (0NC and 1ND). Another difference between both controllers is in their dynamic response, as it is shown in Figure 7.29. It can be seen that in front of an input step of 1 V, the first LQR approach has oscillations of higher amplitude but with lower frequencies, whereas the second approach of cancellation of poles and zeros has higher frequency components but it converges faster to an operating point. Figure 7.29: Comparison of the step response of both LQR and cancellation of p/z approaches. When the observer is implemented together with the LQR controller, the step response is also modified, as the DC gain is affected and the bode diagram slightly increases its peaks due to Modelling and control of the vibration of a rotatory disk p. 101 the dynamics of the Kalman filter. On the other hand, when the low-pass filter is introduced in the cancellation of poles and zeros approach, it is also modified the step response, as the highfrequency components are attenuated. Hence, in Figure 7.30 it is shown the comparison of the final versions of both approaches, where we can see that the second approach has a very fast step response, as the first mode shapes are cancelled and the other ones are attenuated, whereas the first one oscillates more due to the dynamics introduced by the Kalman filter. Figure 7.30: Comparison of the step response of both LQG and cancellation of p/z approaches. The difference between the step response of both approaches is due to the frequency response of each one, compared in Figure 7.31. As presented previously, with the first LQR/LQG approach the amplitude of each peak of the system is attenuated, whereas using the second approach the effect of the first mode shapes is completely canceled, attenuating the remaining ones using a low-pass filter. Figure 7.31: Comparison of the "Next" frequency response using both control approaches in air. p. 102 Report In the case of the controller submerged in water, the step response of the final versions of both approaches is presented in Figure 7.32, where it can be seen that the oscillations are bigger in the case of the LQG, as there are more frequency components as before. However, in the second approach the oscillations are as in the case of the system in air, as the first poles have been cancelled and the last ones are attenuated with the filter. In this particular case, the DC gain of both approaches changes a little bit, so their step response do not converge towards the same value, as it can be seen in the figure. Figure 7.32: Comparison of the step response of LQG and cancellation of p/z approaches in water. Regarding the frequency response, it is shown the comparison of both approaches in Figure 7.33. It can be seen how the different approaches modify differently the frequency response, one attenuating all the peaks of all the natural frequencies and the other one completely cancelling two peaks and attenuating the remaining ones. Figure 7.33: Comparison of the "Next" frequency response using both control approaches in water. Modelling and control of the vibration of a rotatory disk p. 103 As it can be seen, each controller has its own benefits and drawbacks, so in order to choose one to implement in the laboratory it is going to be focused on their performance in a structural way. Therefore, considering that the first modes shapes are the ones that have higher gains and are the ones to be eliminated to improve the performance of the structure, it is better the second approach. Hence, with this second controller it will be possible to fully control the vibration of a specific natural frequency of the disk, rather than attenuating all the vibrations of the disk and not being able to completely cancel any of it, as happens with the first one. In the next section, another control approach is presented to completely reject any type of disturbance entering to the control action. Therefore, the second control approach of the cancellation of poles and zeros is going to be added to that disturbance rejection control to additionally track a "safe" strain reference. p. 110 Report The main drawback of this method is that the part of the controller udthat is in charge of performing the inverse of the plant G−1(s)will have a very high order, corresponding to the same order of the model, which is of order 8 for the case of the system in air and order 12 for the case of the system submerged in water. To avoid having that high order controller, another variation of the same DOB is studied, which consists on using a reduced model for designing the controller udand treat the model mismatch as an internal disturbance wi, which is known as Robust DOB [4]. Therefore, in order to implement the Robust DOB, the plant nominal model G(s)is split in two parts: one corresponding to the first harmonics of the system Gfirst(s), which is a reduced order model used to design the controller ud, and the other part corresponding to the remaining harmonics ∆G(s), which is defined as the model mismatch of the plant ∆G(s) = G(s)−Gfirst(s), introduced in the form of an additive uncertainty, as presented in Figure 8.7. Figure 8.7: Scheme of the Robust DOB. Image from: [4]. By defining this new scheme, the disturbances wconsidered in this case are: the internal disturbance wi= ∆G(s)·u, which represents the model mismatch, and the external disturbance we, considered at the control action. In equation (8.8) both disturbances are presented, together with the computation of this internal disturbance wi. w=we+wi wi= (G(s)−Gfirst(s))u(8.8) In this case, the Robust DOB is computed with the inverse of the reduced plant Gfirst(s), which has the order of the first harmonics to be cancelled (order 4). The weights selected in this case are defined using the filter Q(s), which needs to have a relative degree bigger than the relative degree of the plant model G(s)to be physically implementable (it has to be causal). Due to the specification of reject disturbances at low frequencies, the filter Q(s)is chosen to be a low-pass filter, being the same one as in the case of the implementation of the first DOB approach in water. Modelling and control of the vibration of a rotatory disk p. 111 The Simulink implementation of this control technique is shown in Figure 8.8, where the weight Q2(s)is computed with the reduced version of the plant Q2(s) = Q(s)G−1 first(s), which only includes the first 2 poles of the plant model G(s). It is considered that this plant model G(s) includes the simplified plant Gfirst(s)and the model mismatch ∆G(s)as an additive uncertainty, so the internal disturbance wiis considered into the model, as presented in the previous scheme of Figure 8.7. Additionally, the external disturbance weis considered as the external disturbance that enters in the control action uat the time t= 1 s. Figure 8.8: Simulink implementation of the DOB robust approach using some filters. In this case, it has to be designed the filter Q(s)to accomplish Robust Stability against model misfits, which is checked applying the small gain theorem. Therefore, different block algebra is performed to isolate the uncertainty ∆G(s)from the rest of the blocks, seeking for the equivalent expression of U Wi. This expression corresponds to the equivalent transfer function between the internal disturbance wiand the control action u. Once obtained the U Wiit can be applied the small gain theorem between ∆G(s)and the reduced block, obtaining a bound for the design of the filter Q(s)presented in equation (8.9), which corresponds to the conservative Robust Stability condition.      ∆G(s)−(Hfilter(s)Ccancel(s) + Q(s)G−1 first(s)) 1 + Gfirst(s)Hfilter(s)Ccancel(s)    ∞ <1(8.9) In the case of the disk in air, the Robust Stability condition is not satisfied using the weight of Q(s) = 1 from the previous DOB approach. Hence, a filter that satisfies the inequality of equation (8.9) has to be designed in this case for the Robust DOB in air. We saw that the same low-pass filter Q(s)used in the previous DOB in water fulfilled the Robust Stability inequality for this case in air, so the same filter Q(s)is used, previous presented in equation (8.5). In Figure 8.9 it is shown the performance of this new Robust DOB controller designed for the case of the disk in air, where the udcontrol action is computed using the filter Q(s)and the reduced order plant Gfirst(s), according to Q1(s) = Q(s)and Q2(s) = Q(s)G−1 first(s). On the other hand, if the system is submerged in water, the Robust Stability condition of equation (8.9) has also to be checked considering the reduced system Gfirst(s)and the previous low-pass filter Q(s), obtaining a value of H∞=0.3304. As the infinity norm of that expression is below 1, we can conclude that the system is Robust Stable with that controller and with the uncertainty of the plant ∆G(s), so it can be designed the new Robust DOB using that configuration to cope with the external disturbance and the model mismatch of the plant when the disk is in water. p. 112 Report Figure 8.9: Robust DOB performance in air. In Figure 8.10 it is shown the behavior of the robust DOB in water when tracking a "safe" strain signal and introducing a step input disturbance of 10 V at time t= 3 s. We can see that the Robust DOB takes some time to observe the disturbance and correctly cancels it using the ud control action, as it depends on the filter Q(s). Additionally, it has to be taken into account that the internal disturbance of the model mismatch used in the controller is also present in the simulation, which produces an extra disturbance that is also cancelled by the controller. Figure 8.10: Robust DOB performance in water. Modelling and control of the vibration of a rotatory disk p. 113 The frequency response of this implemented Robust DOB is computed considering that in this case the controller only includes the reduced model of the plant G−1 first(s)with the first 2 poles and 2 zeros of the system. Therefore, when computing the equivalent closed-loop scheme, it is not going to be simplified when multiplied by the plant G(s). This non-simplification effect can be seen in the input-output relation Y/R shown in equation (8.10), which was obtained performing different block algebra. Y R=Q1=Q(s) −−−−−−−−−−−−−−→ Q2(s) = Q(s)G−1 first(s) =G(s)Hfilter(s)Ccancel(s) 1−Q(s) + G(s)Q(s)G−1 first(s) + G(s)Hfilter(s)Ccancel(s)(8.10) Applying this formula, the following frequency response presented in Figure 8.11 is obtained for the Robust DOB. As well as in the case of the standard DOB, it can be seen that even though the filter used Q(s)is not equal to the unity, it still has a small affectation on the frequency response of the controlled closed-loop system. (a) Frequency response in air. (b) Frequency response in water. Figure 8.11: Robust DOB frequency response using different filters Q(s). All in all, with the Robust DOB it has been obtained a more powerful and effective method to achieve robustness against disturbances and model uncertainties, which is perfect for reducing the vibrations of the disk whenever it is in air or water. Moreover, the coupling with the control approach 2 of the cancellation of the first poles and zeros of the plant makes the controlled system converge very fast to a "safe" strain value, neglecting the influence of the first vibrational modes. Modelling and control of the vibration of a rotatory disk p. 115 9 Project Plan The project plan is depicted in the Gantt diagram of Figure 9.1, in which the different tasks performed in this project are labeled in the left and are assigned to the different weeks scheduled for the project. As it can be seen, the project is comprised between February of 2021 and September of 2021, having a total of 28 weeks available. Figure 9.1: Gantt diagram representing the project plan. Modelling and control of the vibration of a rotatory disk p. 117 10 Budget The budget for this project is presented in Table 10.1, in which the different costs of the project are listed. The technical service basically represents the salary of the engineer who is in charge to develop the project, who has an average salary of 15 €/h. Considering the duration of the project of 28 weeks, and assuming an average of worked hours of 28 h/week, the total of worked hours are 784 h, resulting a total amount of budget for only the technical engineer of 11,760 €. Additionally, this project has been supervised by two doctors of the university, so the hours spent by them to guide the project and help with the incoming problems are also taken into account. Other support services such as the office material, the electricity and the internet connection have also been included in the budget for the project. Regarding the office material, 150 €are considered for the whole duration of the project. The electricity consumption is computed considering the total number of worked hours in the computer and the electricity consumed in the lighting and to power all the sensors of the laboratory. In average, it has been obtained an approximated cost of 0.306 €/h, ending up with a total cost of the electricity of 240 €for the whole project. The internet connection is directly computed considering all the months of the project, as the internet is charged monthly independently of the consumption, ending up with a total cost of 210 €. Table 10.1: Total budget for the project. Unit Cost per unit Total cost Technical and support services Technical engineer 784 h 15 €/h 11,760 € Doctors of the university (supervisors) 50 h 40 €/h 2,000 € Electricity consumption 784 h 0.306 €/h 240 € Internet connection 7 months 30 €/month 210 € Research Material ANSYS license software 1 year 4,000 €/year 4,000 € MATLAB license software 1 year 800 €/ year 800 € LabVIEW license software 1 year 3,523 €/ year 3,523 € Research Equipment Amortization of the laptop 7 months 275 €/year 160.4 € Office material - - 150 € Total: 22,843.32 € Industrial profit (5%) +1,142.17 € Total without taxes: 23,985.49 € VAT ("IVA") 21 % +5,036.95 Total budget: 29,022.44 € The software used in this project is contracted by annual subscription, so the full year license has to be paid even if the project only lasts for 7 months. ANSYS multiphysics software has the most expensive license, costing each year 4,000 €/year, while a full LabVIEW license costs 3,523 €/year, and MATLAB costs 800 €/year. The National Instruments equipment and the cost of building the test rig in the laboratory are not considered in the budget for this project, as it is considered that they are already available in the laboratory. p. 118 Report The amortization of the computer is calculated considering a 25% of annual amortization. As the computer was bought by 1,100 €, it results an amortized cost of 275 €/year. However, the project only lasts 7 out of the 12 months of a year, so the total amortization of the computer for the duration of the project is 160.4 €. Finally, after summing all the costs related to the technical and support services, the research material and the research equipment, an industrial profit of the 5% is added to the total budget, which corresponds to an extra amount of 1,142.17 €. Then, the 21 % of the Value Added Tax ("IVA" in Spanish) is added to the resulting quantity, adding an extra 5,036.95 €for the taxes, ending up with the total budget amount for developing this investigation project of 29,022.44 €. Modelling and control of the vibration of a rotatory disk p. 119 11 Environmental impact Regarding the environmental impact of the project, the electrical consumption of the sensors and the computer used for implementing the controller have to be taken into account. In average, a laptop typically uses between 50 and 100 W/hour of electricity when it is being used, depending on the model. Concretely, the Lenovo Legion laptop used for this project has a power consumption of 57 W/hour. Considering that the laptop is used the estimated 784 hours worked for this project, this represents a power consumption of around 44.69 kWh for performing the project, which also includes the time spent in the laboratory performing experiments. Depending on the type of fuel or energy source used to produce this electricity, as well as depending on the type and efficiency of electric power plants, it will be produced more or less amount of carbon dioxide (CO2) to produce this electricity. As an averaged value, it is considered that in September 2021 is produced around 200 g of equivalent CO2per kWh of electricity in Spain, value obtained from the "electricityMap" interactive site. Therefore, to produce the 44.69 kWh needed for the laptop, it is going to be produced around 8.94 kg of CO2. As no natural light arrives in the laboratory, it is needed a constant lighting to illuminate it and work in the test rig. This is considered to be supplied by a 40-Watt bulb equivalent compact fluorescent. Considering that in the laboratory is spent around 5 hours per week, this ends up with a total of light consumption of 5.6 kWh, which is equivalent to produce, using the same CO2emission factors, 1.12 kg of CO2. All in all, the total amount of CO2emitted for lighting the room and powering the computer is estimated to be around 13.880 kg of CO2. To perform numerical experiments in the laboratory, apart from powering all the sensors and the computer used for the controller, it is also needed to introduce water into the tank. The water introduced in the reservoir, which has a maximum capacity of 122 l of water, is not going to be polluted in the experimental process, so it can be safely dumped in the conventional sewage system without the need to treat it. Finally, once the life span of the sensors and the computer used in the test rig ends, they have to be treated as Waste from Electrical and Electronic Equipment (WEEE). The EU has different rules of how to treat these WEEE to address environmental and other issues caused by the growing number of discarded electronics in the EU, which have some materials that are hazardous. Hence, these WEEE have to be collected and treated in order to recycle them, separating the hazardous components and reusing some materials to create new electronic equipment. p. 126 Report [16] H. Schmucker, F. Flemming, and S. Coulson, “Two-way coupled fluid structure interaction simulation of a propeller turbine,” IOP Conference Series: Earth and Environmental Science, vol. 12, p. 012011, 09 2010. [17] L. Learning, “12.3 stress, strain, and elastic modulus.” [Online]. Available: https: //courses.lumenlearning.com/suny-osuniversityphysics/chapter/12-3-stress-strain-an d-elastic-modulus/ [18] CADFEM GmbH, “From ANSYS to System Level Simulation: MOR for ANSYS,” 2008. [Online]. Available: https://www.cadfem.net/fileadmin/user_upload/05-cadfem-in forms/resource-library/From_ANSYS_to_System_Simulation-Model_Reduction_insid e_ANSYS.pdf [19] T. Bechtold, E. Rudnyi and J. Korvink, “mor4ansys: Efficient model order reduction of finite element eletro-thermal mems models.” [Online]. Available: http://modelreduction .com/doc/papers/bechtold05mstkongress.pdf [20] ANSYS Innovation Courses, “Governing equations of modal analysis — lesson 2.” [Online]. Available: https://courses.ansys.com/index.php/courses/modal-analysis/les sons/governing-equations-of-modal-analysis-lesson-2/ [21] MDPI, “Effect of boundary conditions on fluid–structure coupled modal analysis of runners.” [Online]. Available: https://www.mdpi.com/2077-1312/9/4/434/htm [22] ANSYS Inc., “Harmonic analysis in ansys.” [Online]. Available: https://www.mm.bme.h u/~gyebro/files/ans_help_v182/ans_thry/thy_anproc4.html#lspwf3a9mlg [23] ——, “Formulation of Harmonic Analysis in ANSYS,” 2020. [Online]. Available: https://courses.ansys.com/wp-content/uploads/2019/05/3.6.2-Formulation-of-Harm onic-Analysis-New-Template.pdf [24] D. Valentin, A. Presas, E. Egusquiza, and C. Valero, “Influence of the added mass effect and boundary conditions on the dynamic response of submerged and confined structures,” in IOP Conference Series: Earth and Environmental Science, vol. 22, 12 2014, p. 032042. [25] National Instruments, “Measuring Strain with Strain Gages.” [Online]. Available: https://www.ni.com/es-es/innovations/white-papers/07/measuring-strain-with-strai n-gages.html [26] E. Sariyildiz, R. Oboe, and K. Ohnishi, “Disturbance observer-based robust control and its applications: 35th anniversary overview,” IEEE Transactions on Industrial Electronics, vol. 67, no. 3, pp. 2042–2053, 2020. [27] PI Ceramic, “Piezoelectric materials,” February 2021. [Online]. Available: https: //www.piceramic.com/en/products/piezoelectric-materials/ [28] SDTools, “Basics of piezoelectricity,” February 2021. [Online]. Available: https: //www.sdtools.com/help/pz_basics.html#tab%3Apic255FE [29] A. et al., “Negative poisson’s ratio and piezoelectric anisotropy of tetragonal ferroelectric single crystals,” Journal of Applied Physics, 12 2012. [30] OnScale, “Practice: How to calculate piezoelectric material properties from a material datasheet,” February 2021. [Online]. Available: https://onscale.com/blog/practice-how- to-calculate-piezoelectric-material-properties-from-a-material-datasheet/ Modelling and control of the vibration of a rotatory disk p. 127 [31] eFunda Inc., “Piezoelectric transformation between constitutive forms,” February 2021. [Online]. Available: https://www.efunda.com/materials/piezo/piezo_math/transforms .cfm [32] ANSYS Learning Forum, “Rotational velocity,” September 2019. [Online]. Available: https://forum.ansys.com/discussion/9861/rotational-velocity [33] ANSYS inc., “Customer training material, introduction to ansys fluent,” March 2021. [Online]. Available: https://imechanica.org/files/fluent_13.0_lecture05-solver-settings .pdf [34] ANSYS Inc., “Realizable k-epsilon model,” April 2021. [Online]. Available: https: //www.afs.enea.it/project/neptunius/docs/fluent/html/th/node60.htm [35] ——, “ANSYS Command Reference.” [Online]. Available: https://www.mm.bme.hu/~g yebro/files/ans_help_v182/ans_cmd/Hlp_C_CmdTOC.html p. 128 Report Appendices A Numerical modelling of the piezoelectric body A.1 Piezoelectric effect The piezoelectric effect is the property of a material to display electric charge on its surface (producing a change on its polarization) under the application of an external mechanical stress. On the contrary, the converse piezoelectric effect is the production of a mechanical strain due to a change in the polarization of the material. In order to characterize the behavior of a piezoelectric patch, it is needed to define several material properties that define the relationship between the material polarization and its deformation. This relationship can be defined using two different ways: the strain-charge form or the stress-charge form [5]. The strain-charge form is written as: S=sET+dTE D=dT +0rTE(A.1) where S[m m]is the strain, T[Pa]is the stress, E[V m]is the electric field and D[C m2]is the electric displacement field. The material parameters sE[m2 N],d[C N]and rTare the material compliance, coupling properties and relative perimittivity at constant stress. 0is the permittivity of free space, and it has a value of 8.854 ×10−12[F m]. With this notation, the above equations can be written as: Figure A.1: Strain-charge form representation of a piezoelectric material. Image from: [5] On the contrary, stress-charge form relates the physical properties as follows: T=cES−eTE D=eS +0rSE(A.2) The material parameters cE[N m2],e[N V m ]and rScorrespond to the material stiffness, coupling properties and relativity permittivity at constant strain. Modelling and control of the vibration of a rotatory disk p. 129 Within this notation, the above two equations can be written as: Figure A.2: Stress-charge form representation of a piezoelectric material. Image from: [5] A.2 Piezoelectric patch P-876.A12 The piezoelectric patch used in the rotating disc is the Dura Act Piezoelectric Patch Transducer P-876.A12, from PI ceramics. This patch can either work as a sensor or an actuator, as it has the ability to generate and store an electrical charge by means of a laminated structure of piezoceramic, electrodes and a polymer enclosure. The enclosure is used to electrically insulate the patch and make it mechanically robust, as it is shown in Figure A.3. Figure A.3: Composition of the piezoelectric patch from PI ceramic, as well as the polarization direction of each model. Image from: [6] The P-876.A12 model, according to the brochure of the Dura Act Piezoelectric Transducers of PI ceramics (Figure A.4), works with an operating voltage that ranges from -100 to +400 V (thus, making it necessary to use an amplifier to achieve the +400 V), presents a maximum blocking p. 130 Report force of 265 N and is composed by the piezoelectric material PIC255, which is based on modified lead zirconate titanate (PZT) and barium titanate. Figure A.4: Composition of the piezoelectric patch from PI ceramic, as well as the polarization direction of each model. Image from: [6] A.3 Piezoceramic material PIC255 PZT materials can be classified in “soft” and “hard” PZT, depending on the mobility of the dipoles and, hence, on the polarization and depolarization behavior. On the one hand, "soft" PZT materials are used for sensor and actuator applications, as they can be easily polarized even at low field strengths (they present large piezoelectric charge coefficients, moderate permittivities and high coupling factors). Some other applications of soft PZT materials are used for micropositioning, sensors as conventional vibration detectors, ultrasonic transmitters and receivers, or for flow or level measurement, object identification or monitoring and also as sound pickups on musical instruments. The piezoelectric materials classified as soft PZT are the PIC151, PIC255, PIC155, PIC153, or PIC152. On the other hand, "hard" piezoelectric materials can be subjected to high electrical and mechanical stresses, where their properties hardly change under these conditions. Therefore, hard PZT are used in high-power acoustic applications that include ultrasonic cleaning, machining of materials (ultrasonic welding, bonding, drilling,..), ultrasonic processors, medical sector, and sonar technology. Some of these materials that are categorized as hard PZT are PIC181, PIC184, PIC144, PIC241, PIC300 [27]. Modelling and control of the vibration of a rotatory disk p. 131 As stated in the brochure information of Figure A.4, the material used in the DuraAct patch is the PIC255, which is defined as a soft PZT. The piezoelectric constants provided by the material datasheet of PI Ceramics are presented in Figure A.5. Figure A.5: Specific parameters of the standard materials, from PI Ceramic. Image from: [7] p. 132 Report The constants given by the manufacturer represent the following physical properties [8]: •Curie Temperature Tc[◦C] The Curie Temperature is defined as the temperature at which the crystal structure of the material changes from piezoelectric structure (non-symmetrical) to non-piezoelectric structure (symmetrical), loosing all its piezoelectric properties. •Permittivity constant  The permittivity, also called dielectric constant , is the dielectric displacement per unit of electric field. Depending on the superscript of this constant it means that it is measured at constant stress Tor at constant strain S. The subscripts ij correspond to ithe direction of the dielectric displacement, and jthe direction of the electric field. The dielectric constant, as the case of the datasheet, is usually given as a ratio of the permittivity of the material to the permittivity of free space 0= 8.854 ×10−12F/m. The values given by the manufacturer for PIC255 are T 11 = 1750 ·0, which is the permittivity for dielectric displacement and electric field in direction 1 (perpendicular to direction in which ceramic element is polarized), under constant stress, and T 33 = 1850 ·0, which is the same but in direction 3 (parallel to direction in which ceramic element is polarized). •Dielectric Loss Factor tan δ The dielectric loss factor, also known as dissipation factor, is defined as the tangent of the dielectric loss angle tan δ, which determines the ratio of effective conductance to effective susceptance in parallel circuit, measured by using an impedance bridge. The values for tan δ are usually determined at 1 kHz. •Electro-mechanical coupling factor k This factor is defined as the ratio of the mechanical energy accumulated in response to an electrical input or vice versa. It is an indicator of the effectiveness with which a piezoelectric material converts electrical into mechanical energy, or converts mechanical into electrical energy. The first subscript of k denotes the direction along which the electrodes are applied, whereas the second subscript denotes the direction along which the mechanical energy is applied. •Piezoelectric charge coefficient d It represents the polarization generated per unit of mechanical stress T [C/N] applied to a piezoelectric material or, alternatively, is the mechanical strain S experienced by a piezoelectric material per unit of electric field applied [m/V]. The first subscript indicates the direction of polarization generated in the material when the electric field, E, is zero. The second subscript is the direction of the applied stress. Alternatively, the first subscript can be the direction of the applied field strength and the second one the induced strain. •Piezoelectric voltage coefficient g It is the ratio of the electric field generated by a piezoelectric material per unit of mechanical stress applied [Vm/N] or, alternatively, is the mechanical strain experienced by Modelling and control of the vibration of a rotatory disk p. 133 a piezoelectric material per unit of electric displacement applied. The first subscript indicated the direction of the electric field generated in the material (or the direction of the applied electric displacement), whereas the second subscript is the direction of the applied stress (or the induced strain). This factor is important for assessing a material’s suitability for sensing applications, as it depends on the strength of the induced electric field produced by a piezoelectric material in response to an applied physical stress. •Acoustic-mechanical frequency coefficient N The frequency coefficient Nis the product of the resonance frequency fsand the linear dimension of the ceramic element governing the resonance. Depending on the geometry of the patch it can be defined 3 different constants: the Radial Mode Disc resonance frequency constant Np=fs·Dφ(where Dφis the diameter of the ceramic element), the Thickness Mode Disc frequency constant NT=fs·h(where his the thickness hof the ceramic element), and the Longitudinal Mode frequency constant NL=fs·l(where lis the length of the element). •Elastic compliance coefficient s Is the strain produced in a piezoelectric material per unit of stress applied and, for the 11 and 33 directions, is the reciprocal of the modulus of elasticity (Young’s modulus Y). We can distinguish between sD, which is the compliance under a constant electric displacement D, and sE, which is the compliance under a constant electric field E. The first subscript indicates the direction of strain and the second one the direction of stress. Notice that, as the piezoelectric is an anisotropic material, all the piezoelectric constants provided by the manufacturer present two subscripts that relate the direction of the electrical field associated with the voltage applied, or the charge produced, with the direction of mechanical stress or strain. Superscripts indicate a constant mechanical or electrical boundary condition. Normally, the direction of the positive polarization is made to coincide with the Z-axis of the rectangular system (subscript 3), where directions X and Y are represented by the subscript 1 and 2, respectively. Shear about these X, Y, Z axis are represented by the subscripts 4, 5, 6, respectively, as schematized in Figure A.6 [8]. Figure A.6: Directions of forces affecting a piezoelectric element. Image from: [8] p. 134 Report A.4 Computation of the constants demanded by ANSYS Using all the different constants given by the manufacturer for the PIC255 material, it is needed to transform all these data into a format which can be inserted into the material database of ANSYS, as well as in the properties of the piezoelectric body to simulate the piezoelectric behavior. As stated before, the elastic compliance coefficient sis the reciprocal of the Young’s Modulus, so the modulus Y11,Y22 and Y33 can be computed from these values: Y11 =Y22 =1 sE 11 =1 16.1×10−12 m2/N = 62.112 GPa (A.3) Y33 =1 sE 33 =1 20.7×10−12 m2/N = 48.309 GPa (A.4) Knowing the values of sE 11,d31,T 33,kp,sE 12 we can obtain the exact value of the transversal Poisson’s ratio ν12 [28]: sE 12 =−sE 11 + 2 d2 31 k2 pT 33 =−5.808 ×10−12 m2/N (A.5) ν12 =−Y11sE 12 = 0.3608 (A.6) According to the material datasheet of Figure A.5, all the piezoelectric materials have a Poisson ratio around 0.34, so we can check that the computed value is around there. Then, the shear stress value is computed as: G12 =Y11 2 (1 + ν12)= 22.82 GPa (A.7) sE 55 =d2 15 T 11k2 15 = 4.482 ×10−11 m2/N (A.8) G13 =G23 =1 sE 55 = 22.31 GPa (A.9) According to [29] the Poisson’s ratio of a material is related with the piezoelectric anisotropy following the relation ν13 =−sE 13/sE 33. However, the author states that in practice it can be estimated following the semi-empirical relation of: d33 d31 =−1 ν13 (A.10) Obtaining that the longitudinal Poisson’s ratio (relates the polarization direction of the patch with its perpendiculars) is: ν13 =ν23 =−d31 d33 =180 400 = 0.45 (A.11) Once obtained all the piezoelectric constants needed, we can build the elastic compliance matrix Modelling and control of the vibration of a rotatory disk p. 135 [sE]for the orthotropic material using the following expression [30]: [sE] =           1 Y11 −ν21 Y22 −ν31 Y33 0 0 0 −ν12 Y11 1 Y22 −ν32 Y33 0 0 0 −ν13 Y11 −ν23 Y22 1 Y33 0 0 0 0001 G23 0 0 0 0 0 0 1 G13 0 0 0 0 0 0 1 G12           (A.12) Recalling that the number 3 is referred to the poling direction Z and 1,2 with the directions X,Y respectively we can define the following equivalences due to symmetry: ν21 Y22 =ν12 Y11 ,ν31 Y33 =ν13 Y11 ,ν32 Y33 =ν23 Y22 (A.13) Then, the compliance matrix ends up being: [sE] =           1 Y11 −ν12 Y11 −ν13 Y11 0 0 0 −ν12 Y11 1 Y22 −ν23 Y22 0 0 0 −ν13 Y11 −ν23 Y22 1 Y33 0 0 0 0001 G23 0 0 0 0 0 0 1 G13 0 0 0 0 0 0 1 G12           (A.14) With this matrix we can then change from the strain-charge (A.1) to the stress-charge form (A.2), which is the form that the PiezoAndMems extension of ANSYS uses to characterize the piezoelectric effect. This transformation [31] is done computing the inverse of the previous matrix [cE] = [sE]−1, obtaining: [cE] =          1.1861 0.7296 0.6705 0 0 0 0.7296 1.1861 0.6705 0 0 0 0.6705 0.6705 0.9524 0 0 0 0 0 0 0.2231 0 0 00000.2231 0 000000.2282          ×1011 Pa (A.15) Finally, together with the dcoefficient matrix, it can be computed the [e]matrix by performing [e]=[d]·[cE], obtaining the following non-zero coefficients demanded by ANSYS: e31 =e32 =−7.6628N/Vm (A.16) e33 = 13.9597 N/V m (A.17) e15 =e24 = 12.2716 N/V m (A.18) The units used in ANSYS for expressing all these constants are the SI base units, which are 7 base quantities (s, m, kg, A, mol, cd) from which all other SI units can be derived. Concretely, the coefficients of the [e]matrix are expressed in [A·sec m2]≡[N V·m], and the permittivitiy is expressed in [A2·sec4 kg·m3]≡[C V·m]. p. 142 Report the mesh used for the ANSYS Mechanical transient simulation has a total of 29566 nodes and 12534 mesh elements. In Figure C.1 it can be seen the mesh used, where the zones around the piezoelectric element has a bigger number of elements than the rest of the disk. Figure C.1: Section of the Mesh of the system for the Transient Structural simulation. C.3 Analysis settings In the Analysis Settings of the Transient Structural simulation it is specified several information regarding the time step used in the simulation and different information about the data management. Concretely, it can be defined the duration of the simulation and the different steps that ANSYS has to perform to solve the model, which can be used to increase the accuracy of the dynamic response of the disk. In the case of the 2-way FSI analysis, at each time-step the System Coupling will send back and forth data to the ANSYS Mechanical simulation, so it’s important to define the desired simulation duration in "Step End Time" (corresponding to the total number of seconds to last the simulation), and also to set the "Auto Time Stepping" to to Off, defining as 1 the total number of substeps for the simulation, which will produce hat the the System Coupling exchanges data every 1 time step. C.4 Rotational boundary conditions The disk in the test rig has a motor that spins the disk at a rotational speed up to 3000 rpm, so the inertia effects of the rotating parts must be consistently represented in order to accurately simulate the system. In order to define a rotational velocity in ANSYS Mechanical, it can be done using APDL by means of the OMEGA or CMOMEGA commands. The OMEGA command is used to specify the rotational velocity of the structure about global cartesian axes, used when the whole model is rotating, whereas CMOMEGA is used to specify the rotational velocity of an element component about a user-defined rotational axis, being able to define an stationary part and/or several rotating parts having different rotational velocities. It is important to mention that the rotation in the transient structural analysis of ANSYS it is not applied directly by actually spinning the model, but just as a radial acceleration load to each element, which depends on the radius of the element from the axis and the rotational velocity [32]. This fact is one of the problems that appear when performing the 2-way FSI simulation, Modelling and control of the vibration of a rotatory disk p. 143 because it is needed to spin the model to have a direct effect on the fluid flow analysis, so the radial acceleration load that ANSYS applies to each element when defining a rotational velocity is not enough for the case. For that reason, in the 2-way FSI simulation it is needed to define directly a rotational velocity to all the solids in order to spin them and have a direct affectation of the rotation to the fluid. This boundary condition can be applied as a nodal velocity with the APDL commands IC and ICROTATE, used to generate equivalent nodal translational velocities for lines and solid elements that rotate about an axis. However, defining a rotational velocity in the transient simulation, and waiting for the coupled 2-way FSI system to converge with the rotation of the fluid, it is a very slow and prone to errors process, so it has been decided to apply the equivalent rotation to the solids using the aforementioned CMOMEGA command. By doing so, the effect of the rotation has to be also introduced in the Fluent simulation of the coupled analysis, as it is not taken into account in the transient simulation. This way, the convergence of the system is faster and it does not yield to unexpected errors when synchronizing both simulations in the 2-way FSI. To sum up, the rotation in ANSYS Mechanical is defined by the APDL command of CMOMEGA, which is equivalent to use the Rotational velocity boundary condition, as it only accounts for the structural effects of a part spinning at a constant rate. In Figure C.2 it can be seen an overview of the rotation applied to the shaft, as well as the boundary condition applied to account for the bearing effect. Figure C.2: Rotational velocity definition in ANSYS Workbench, together with the cylindrical support. p. 144 Report Together with the rotational velocity boundary condition, it has also to be defined the constraints of the bearings, which do not allow the shaft to move neither longitudinally nor transversely. This boundary condition corresponds to set both faces of the bearings to be cylindrical supports, allowing the solid to freely rotate tangentially but constraining the radial and axial movement of it. C.5 Piezoelectric boundary conditions One challenging part of this simulation is the integration of the piezoelectric element, which bends and applies an stress when a voltage is applied to it. However, this electromechanical coupling is handled by the PiezoAndMems extension, in which it is defined the the active layer of the piezoelectric patch as a "Piezoelectric Body". In this definition, it has to be inserted all the piezoelectric properties in the suitable ANSYS format, which is computed in Appendix A.4 from the data given by the piezoelectric manufacturer, and is presented in Table C.1. Moreover, the polarization axis has to be defined, which in the case of the patch used in the disk corresponds to the axis perpendicular to the surface of the transducer, the Z axis. Table C.1: Properties of the "Piezoelectric Body" defined in ANSYS. Permittivity of free space 0= 8.854 ×10−12 Piezoelectric stress constant PIEZ e31 =−7.6628 N/Vm e33 = 13.9597 N/Vm e15 = 12.2716 Dielectric permittivity DPER 11 = 1750 33 = 1850 Once defined the piezoelectric body, it has to be applied the corresponding boundary conditions of the voltage to each electrode of the active layer. Knowing the composition of the P-876.A12 patch presented in Appendix A.2, it is applied a positive voltage to the upper face of the patch and a negative voltage to the lower part, setting the sign of the polarization and creating the electrical boundary condition for the piezoelectric patch that can be seen in Figure C.3. (a) "Piezoelectric Body" assigned to the patch. (b) Voltage boundary condition. Figure C.3: Details of the piezoelectric defined in ANSYS using the PiezoAndMems extension. Modelling and control of the vibration of a rotatory disk p. 145 C.6 Fluid-Solid Interface (only in the 2-way FSI simulation) This boundary condition of the Transient Structural simulation is only used in the case of the 2-way FSI simulation, as it defines the surfaces of the solid that are in contact with the fluid, which are the ones that will send displacement data to the fluid analysis and will also receive force data back from the fluid analysis. This surface, known as "Fluid-Solid Interface", includes all the outer walls of the rotating structure that are in contact with the surrounding water of the tank, which include to the outer faces of the passive layer of the piezo, the walls of the disk, and the part of the shaft that is submerged in the water. The faces that are not part of this "Fluid- Solid Interface" are the ones that are not in direct contact with the water, which are the face between the disk and the piezoelectric and the face between the shaft and the disk. Once defined this interface, the Transient Structural part of the 2-way FSI is fully defined, as presented in Figure C.4. Figure C.4: Boundary conditions of the Transient Structural Simulation. Modelling and control of the vibration of a rotatory disk p. 147 D Fluid Flow (Fluent) detailed simulation The Fluid Flow (Fluent) system is able to perform fluid-flow analysis of incompressible and compressible fluid flows, as well as being able to simulate laminar and turbulent flows, using ANSYS Fluent, which provides comprehensive modeling capabilities together with the capacity of solving conservation equations for mass and momentum. This analysis is configured in ANSYS Fluent from an imported computational mesh of the geometry, which is used by Fluent to define pertinent mathematical models, select materials, define boundary conditions, and specify solution controls that best represent the problem to be solved. Finally, after Fluent solves the mathematical equations, the results of the simulation are displayed in CFD-Post for further analysis (for example, contours, vectors, and so on). D.1 Geometry and Mesh Analogously as in the Transient Structural simulation of Appendix C, first of all it is needed to suppress the solids of the geometry that are not involved in this analysis. Concretely, all the solids unless the tank of water "tank" have to be suppressed, as only the fluid is considered in this analysis. The properties of the water inside the tank are not needed to be defined as a material for the solid, as the type of fluid will be selected in the following Fluent interface. Regarding the mesh of the tank of water, a selective body meshing technique is used with a patch conforming method to set up a tetrahedral mesh inside the tank of water. Tetrahedral mesh elements are used to generate mesh on the cylindrical domain with an element size of 0.01 m. Smoother mesh elements are generated around the shaft and disk surface to capture the boundary layer flow. Then, an inflation around the interior faces is used, defined with the Total Thickness inflation option. The maximum thickness defined is set to 0.01 m, with a growth rate of 1.2 and the total number of layers of 5. Doing so, the resulting mesh is the one shown in Figure D.1, where it can be clearly seen that the tank does not include neither the disk nor the piezoelectric element, and it can also be seen the inflation of the mesh near the disk to account for the boundary layer flow. This mesh has a total of 183,210 nodes and 97,1335 mesh elements. Figure D.1: Section of the mesh used in fluent analysis. p. 148 Report For defining the boundary conditions in the ANSYS Fluent application, it is also needed to create some named selections in the meshing application. This named selections correspond to assign the top face of the tank with a name "opening", whereas the rest of the outer faces of the cylinder are named "walls". Inside the tank of water, the walls that are in contact with the disk are also named as "contact_disk", and are the ones that will be receiving and sending back information to the Transient Structural application. The mesh is then automatically imported into the Fluent interface by updating the corresponding setup cell of the Fluid Flow (Fluent) analysis. The Fluent interface is presented in Figure D.2, in which the imported geometry is displayed in the right hand side of the window, with the possibility to set the front faces transparent to see the interior faces, and in the left hand side of the windows is placed the Outline View, where all the properties for the Fluent analysis have to be defined. Figure D.2: Fluid Flow (Fluent) interface. In this Fluent interface two important blocks have to be configured, which are marked in the Outline View: the Fluent Setup, where the material and the boundary conditions are set, and the Fluent Solution, in which specific solution parameters are specified and the solution is initialized. D.2 Fluent Setup General First of all, various generic problem settings are defined in the General tab of the Fluent Setup, such as those related to the mesh or the solver. In the "Mesh" section of the Task Page, which can be seen in the previous Figure D.2 between the Outline View and the Model, it can be displayed and checked some controls related to the mesh settings such as verifying the validity of the mesh with some volume statistics, mesh topology and periodic boundary information, as well as displaying the quality of the mesh. All those mesh parameters play a significant role in the accuracy and stability of the numerical computation of the system, so it’s important to check them. Modelling and control of the vibration of a rotatory disk p. 149 Moreover, some controls relating the "Solver" settings can also be defined in the same General tab. In the Solver "Type" two different solution methods are available for computing a solution for the model: the Pressure-Based and the Density-Based. The Pressure-Based analysis enables the pressure-based Navier-Stokes solution algorithm, extracting the pressure field by solving a pressure equation which is obtained by manipulating continuity and momentum equations. On the other hand, Density-Based analysis enables the density-based Navier-Stokes coupled solution algorithm, where the continuity equation is used to obtain the density field. In the literature says that historically, the pressure-based solver was developed for low-speed incompressible flows, while the density-based approach was mainly used for high-speed compressible flows, being more accurate for transient cases with combustion, hypersonic flows or travelling shocks. Therefore, as the pressure-based solver is applicable for flow regimes from low speed incompressible to high-speed compressible flows, requires less memory (storage) and allows flexibility in the solution procedure, it has been selected as the solver type for the Fluent analysis [33]. Additionally, the "Velocity Formulation" specifies the formulation to be used in the velocity’s calculation, and can either be Absolute or Relative. In some manuals is recommended to use the velocity formulation that will result in most of the flow domain having the smallest velocities in that frame, thereby reducing the numerical diffusion in the solution and leading to a more accurate solution. Absolute velocity formulation is normally preferred in applications where the flow in most of the domain is not moving, whereas the Relative velocity formulation is appropriate when most of the fluid in the domain is moving, as in the case of a large impeller in a mixing tank. However, when a coupled solution algorithm is used, the relative velocity formulation is not available, and only the absolute formulation can be used. Regarding the "Time" dependence solution, as it will be synchronously solved with the Transient Structural, it is selected to be performed a transient analysis in Fluent. Finally, the "Gravity" is enabled and specified its value in the corresponding axis of the vertical orientation of the object, corresponding to the Z axis in the "Gravitational Acceleration" tab. Models Different type of turbulence models can be defined for the Fluid analysis, such as Multiphase, Energy, Viscous, Radiation, Heat Exchanger, Species, Discrete Phase, Acoustics or Structure. For the fluid analysis of the water inside the tank, "Viscous Model" is selected, in which three different k−models can be specified to calculate the turbulent flow (Standard, RNG or Realizable). In this case, the Realizable k−model has been selected, which computes the dissipation rate using a derived transport equation from the exact equation for the transport of the meansquare vorticity fluctuation, using a variable turbulent viscosity Cµ[34]. This difference in the computation of the dissipation rate produces that Realizable k−model has a superior ability to provide superior performance for flows involving rotation. Additionally, a "Scalable Wall Function" model is selected to avoid the deterioration of standard wall functions under a small grid refinement, producing consistent results for grids of arbitrary refinement. p. 150 Report Materials By default, Fluent has defined the air as the Fluid material for the fluid simulation, so it is needed to change it to water liquid for the case study. Checking into the Fluent Database Materials, it can be selected the "water-liquid(h2o<l>)" material, which is the one representing the fluid around the disk of water, and it has already defined all the properties of water needed for the simulation. Cell Zone Conditions Once defined the material, in the Cell Zone Conditions tab it is selected that the whole volume of the tank is made of water. Moreover, some configuration of the fluid cell zone can be defined by entering in the Fluid dialog box, which opens the window shown in Figure D.3. As it can be seen in the Figure, it is activated the "Mesh Motion" option, which enables the sliding mesh model for the cell zone. On the configuration of this mesh motion it can be specified the translational and/or rotational velocity of each moving fluid zone. As in the real configuration the tank of water is static and the water is moved due to the rotation of the disk, the approach done in the simulation is not to move the fluid itself but move the walls that are in contact with the disk, which will produce a movement to the fluid as a consequence. Therefore, in the cell zone condition tab is set a rotational and translational velocity of 0. Figure D.3: Fluid dialog box of the "tank" Cell Zone Condition. Boundary Conditions The Boundary Conditions task page allows you to set the type of a boundary and display other dialog boxes to set the boundary condition parameters for each boundary. In Figure D.4 there’s the summary of the Boundary Conditions defined in each wall of the tank of water, which are separated by the different type of boundaries, corresponding to the internal boundary, the outlet Modelling and control of the vibration of a rotatory disk p. 151 and the wall. Moreover, next to the zone name it appears an identifier (id), which is just used by Fluent for informational purposes, so it cannot be edited. Figure D.4: Summary of the boundary conditions applied to each wall of the tank of water. The "internal" boundaries are defined automatically by the volume occupied by the tank of water, and the "outlet" is assigned to the top surface of the cylinder, named previously as "opening". This outlet is defined with a "pressure-outlet" boundary condition, with a relative pressure of 0 Pa, in the direction normal to boundary. To check that the absolute operating pressure is correctly set to 1 atm, it can be done in the Operating Conditions tab, where it can be seen that the atmospheric pressure has a value in the simulation of 101325 Pa. The "contact_disk" walls, which are the faces in contact with the disk and the shaft, are set to type "wall" and are defined as a "Moving Wall" with a relative rotational speed of 6.28 rad/s, which is the same value as the rotational speed introduced in the disk in the Transient Structural simulation. In Figure D.5 it can be seen the details of this "wall" boundary condition for the case of the "contact_disk" zone. Using this approach, the fluid inside the tank of water will rotate accordingly to the rotation of these walls, analogously as the real system will do. Figure D.5: Details of the boundary condition of the walls that are in contact with the disk. The internal faces of the tank of water, named "outer_walls", are also defined as a boundary condition of type "wall", but in this case they are defined as a "Stationary Wall", with a "No Slip" shear condition.