The role of bond graph modeling and simulation in mechatronics systems

Mechatronics 12 (2002) 1271–1295

The role of bond graph modeling and simulation in mechatronics systems An integrated software tool: CAMP-G, MATLAB–SIMULINK Jose J. Granda Department of Mechanical Engineering, California State University, Sacramento, CA 95819, USA Fachhoschule Bonn-Rhein-Sieg, University of Applied Sciences. Sankt Augustin, Germany

Abstract One of the main and most challenging steps in the design and analysis of a mechatronics system is to generate a computer model. This paper explores the fundamental theory, the methodology and the process from conceptual ideas to practical realization. Using a multienergetic approach that allows the modeling of interdisciplinary models, it explores the theory and method to automate the process of the generation of the differential equations and how to automate the derivation of transfer functions. The approach is discussed for linear and nonlinear systems. The generation of a computer model takes new dimensions when that model contains mixed energy domains such as electromechanical, electrohydraulic, thermo-fluid, and electronic control systems all together. These are typical of mechatronics applications. This paper explores the bond graph technique as a modeling tool to generate state space models or non-linear models together with software tools. CAMP-G (Computer Aided Modeling Program with Graphical input) has been developed in order to generate computer models automatically and have them integrated with MATLAB–SIMULINK as simulation tools. Several aspects of mechatronics systems design have been investigated in order to focus on which areas the bond graph modeling technique can help engineers in the process of creating mechatronics systems from scratch. Towards this end, the paper deals with computer-generated models of sensors, actuators, and multidisciplinary complex physical systems. Ó 2002 Published by Elsevier Science Ltd.

E-mail address: [email protected] (J.J. Granda). URL: http://gaia.csus.edu/
grandajj. 0957-4158/02/$ - see front matter Ó 2002 Published by Elsevier Science Ltd. PII: S 0 9 5 7 - 4 1 5 8 ( 0 2 ) 0 0 0 2 9 - 6

1272

J.J. Granda / Mechatronics 12 (2002) 1271–1295

1. Introduction Currently the subject of mechatronic system design has taken a new dimension since many researchers have turned their attention to emerging technologies that allow engineers to understand multidisciplinary systems. The ASME for example, in their publication [7], gives a clear idea of the need to explore new techniques that are to be at the center of mechatronics. One of these is the bond graph modeling technique. Experts in mechatronics interviewed in [7] report why bond graph modeling will play a major role in the design of mechatronics systems in the years to come. Currently, there are over 400 research publication and many books about the bond graph method as evidenced by the existence of the bond graph compendium of literature [8]. For the reader not familiar with bond graph methodology, a good reference is a textbook such as [2]. Bond graph modeling is used in industry, for example the auto industry where many applications contain electromechanical devices, sensors, actuators, and microprocessors driven systems. Many other industries such aerospace, hydraulics, power, chemical and lately the biomedical have the need for methods and tools to design systems and processes across the disciplines. On this point, this paper presents a multidisciplinary approach for these solutions. In the process of creating a physical device from scratch, one of the major tasks is the creation of the physical model, that of its control system, and that of its sensors and actuators in order to allow software tools to be integrated in the design process. This paper explores an integrated approach between theoretical principles of bond graph modeling and software. This includes a graphical software tool to create computer models automatically called Computer Aided Modeling Program with Graphical input (CAMP-G). This tool implements the theoretical principles outlined here and then the models are transformed into MATLAB and the graphical nonlinear simulation tool SIMULINK. Using a bond graph model and then CAMP-G, one can generate first order state space differential equations, symbolic transfer functions and SIMULINK S functions in source code form such as MATLAB .m or .s files, which SIMULINK can process. The fundamentals of such approach are discussed in [9]. The idea here is to generate automatically the differential equations from the bond graph model. This means the generation of the system A, B, C, D matrices to produce the state space representation, it means also the computer generation of transfer functions and output equations in symbolic form. For non-linear systems this implies a set of equations where the individual non-linear constitutive relations can be accomplished as well as any conditional switches, nonlinear time dependent effects and of course a form suitable for logical execution of all these by MATLAB and the non-linear S functions of SIMULINK.

2. Piezoelectric sensors Piezoelectric sensors are used in many engineering applications to measure physical systems variables. The device such as the one shown in Fig. 1 consists of an oscillating mass with damping and stiffness effects enclosed in a housing. The force

J.J. Granda / Mechatronics 12 (2002) 1271–1295

1273

Fig. 1. Piezoelectric sensor.

acting on the mass is transmitted to a piezoelectric material in the form of a wafer attached to the oscillating mass. When such material is loaded, produces an electric charge associated with the mechanical deformation. Since such a charge is small, it requires an operational amplifier in order to be measured and find a relation to the acceleration, which the devise is trying to measure. There are essentially three parts to the sensor. One is the mechanical part, the second is the piezoelectric transformation and the last one is the electronics necessary to produce the amplification. The oscillating mass is a second order mass, damper, stiffness system. Let us consider now the models of each of these sections with the aim at putting them together as a multidisciplinary mechatronics model. 2.1. Operational amplifier model From the bond graph representation shown above, one can make the equivalence of the equations that control the operation of the amplifier obtained from physical principles and those obtained using the bond graph model notation. Basic operation of the operational amplifier requires that the input current be zero. In order to achieve this, the operational amplifier has to have very high input impedance. Representing the voltage at the input node, U2 (Fig. 2a) is a small voltage and in an ideal case the whole current iq flows through the resistor and the capacitor so that the current ia is zero. The current iq which originates in the piezoelectric transformation, can be expressed as: iq ¼ i2 . It is equal to the sum of the currents, which flow through the resistor and the capacitor iq ¼ 

Ua dUa  C R dt

ð1Þ

which leads to the differential equation of the electrical circuit in terms of the output voltage RC

dUa Ua dq þ ¼ R dt dt R

ð2Þ

1274

J.J. Granda / Mechatronics 12 (2002) 1271–1295

Fig. 2. Operational amplifier models: (a) operational amplifier, (b) operational amplifier bond graph model.

Fig. 2b represents the bond graph model of the amplifier obtained in a topological manner. The zero junctions are the points of common potential and they have been laid out to resemble the physical layout of the circuit. These zero junctions are points of distinct voltages. The amplifier gain is represented by a controlled source whose input is only the voltage signal U2 and its voltage amplification is U2 A. Since Ua is a high voltage compared to U2 , it follows that the current flows as indicated through the R- and C-elements. The voltage across the R- and C-elements is the difference between Ua and U2 and the 1 junction represents such. Both the R- and C-elements are attached to the 0 junction to indicate the common potential across them. This also indicates that the current i2 is the sum of the currents through C and R and thus ic þ ir ¼ i2 . The arrows in the 0 and 1 junctions indicate how the currents or the voltages sum or subtract. Fig. 3 indicates an augmented and simplified bond graph where the variables have been generalized as R, C, and SE elements referenced by their bond number and the type of element they represent. Such notation allows the description of a physical model by a bond graph into computer software like CAMPG in such generalized form so that systems in different energy domains can be processed all together. The e variables are voltages and the f variables are currents in this case.

Fig. 3. Operational amplifier simplified bond graph model.

J.J. Granda / Mechatronics 12 (2002) 1271–1295

1275

Using the bond graph notation above, junctions are numbered and the physical parameters are distinguished by the bond number to which they are attached. Using the 0 and 1 connections together with the individual element constitutive relations, the differential equation can be derived in terms of the capacitor charge q as a state variable. The variable chosen here is the charge as a state because in fact it is the current i2 that gets integrated. The corresponding differential equation dq12 =dt is written as dq12 and the entire equation in a format compatible with compilers or interpreter software. The equations below follow the mathematical form known as the Cauchy form. The voltage output equation is e11 or e12 : dQ12 ¼ SF6  Q12=C12=R13 e11 ¼ Q12=C12

ð3Þ

Here one can clearly recognize that f7 ¼ f6 , which in fact is iq ¼ i2 in the schematic of Fig. 3. Also f11 ¼ f7 ¼ f6 ¼ iq and f12 ¼ f11  f13 , f6 ¼ iq , f12 ¼ ic and f13 ¼ ir . Now it can be demonstrated that these equations have a direct relation to those used with the conventional notation. f12 ¼ f11  f13 is ic ¼ iq  ir and thus The following equation is obtained for the current through the capacitor f12 ¼  

dUa C dt

dUa Ua C ¼ iq  dt R

ð4Þ ð5Þ

and thus yields the differential equation below: RC

dUa Ua dq þ ¼ R dt dt R

ð6Þ

This clearly demonstrates that this is the same as Eq. (2) obtained by conventional methods based on the summation of currents at the nodes following Kickoff’s law. Hoffman [6] obtained an equation like the one shown as Eq. (6) applying basic principles. This verifies the validity and complete equivalence of the bond graph representation of the operational amplifier. One should notice that the choice of state variable q in this case (charge) is automatically made to be the result of the physical variable that gets integrated (current I). One can take the voltage Us as state but a lot of manipulation of the equations is necessary such as shown in [6]. In presenting this approach this complexity is simplified and the voltages (e), effort variable, will be calculated. The objective then would be that these would be computer generated in symbolic form. 2.2. Piezoelectric transformation The piezoelectric transformation is a direct relation between the displacement of the mass and the charge generated in the capacitor so that it would yield: q ¼ Kq y

ð7Þ

1276

J.J. Granda / Mechatronics 12 (2002) 1271–1295

Fig. 4. Piezoelectric transformation.

and thus dq dy ¼ Kq dt dt

ð8Þ

which represents the piezoelectric transformation between the velocity of the oscillating mass and the current, which is the input to the operational amplifier. This equation in bond graph notation is represented by a transformer element, which transforms the input force into voltage proportional to the acceleration. This can be represented as shown in Fig. 4. 2.3. Mechanical model The mechanical model consists of the housing and a mechanical oscillator with stiffness and damping. The physical system and bond graph are shown in Fig. 5. Here the velocities y_ and u_ are absolute velocities with respect to an inertial frame. Considering that the piezoelectric effect is produced by the relative displacement of the mass on the piezoelectric material then it makes sense to express the bond graph and indeed the equations of motion of the mass in terms of the relative motion with respect to the housing similarly it is attached to the source of €u. This means, a simplification of the model is useful, without considering the mass of the housing, but considering that the housing transmits acceleration to the mass. The bond graph shown in Fig. 6 represents this simplification. The physical parameters and expressions for the forces have been superimposed on the bond graph for clarity, but they are not required as the bond graph uses the effort variables and the bond numbers to express the generalized forces and generalized velocities. Generalized forces (efforts e), generalized velocities (flows f).

Fig. 5. Sensor mechanical section and equivalent bond graph model.

J.J. Granda / Mechatronics 12 (2002) 1271–1295

1277

Fig. 6. Bond graph in terms of relative motion.

The forces superimposed on the bond graph of Fig. 6b help to demonstrate the equivalence of the differential equations that one can obtain using the free body diagram of the mass, damper spring system. The 1 junction represents the summation of forces in the same way as using Newton’s law: m€yr þ b_yr þ kyr ¼ m€ u ð9Þ Now it is necessary to assemble the piezoelectric sensor as a complete system. 2.4. Electromechanical integrated sensor model and computer simulation All three sections of the sensor explained above, can now be put together as a single multidisciplinary model. This paper is oriented to show the generation of models and the simulation of multidisciplinary systems known as mechatronics system. To this end, one can generate the transfer function between the output voltage and the input acceleration, variables that are of interest in this case. Now we proceed to transform the bond graph model into a computer model. To achieve this, the bond graph is entered in CAMP-G in graphical form using the icons on the left as shown in Fig. 7. Software has been developed to interface CAMP-G with MATLAB, one can obtain computer-generated .m files containing the system equations, the state space matrices and the desired transfer functions. The method employed here is different than the one employed in [6]. In such reference, the author derived several individual transfer functions for each of the three sections of the sensor, the mechanical, the piezoelectric transformation and the electrical part. Then, these were combined later to obtain the voltage/acceleration transfer function. In contrast, here, the model is a multienergy, multidisciplinary model which combines all three sections and treats the model as a whole with one single generalized notation. The effort variables of the mechanical section represent the forces, the effort variables of the piezoelectric transformation represent the relation between the forces, which the sensor is subjected and the voltage produced as a result of the piezoelectric effect. These variables in the electrical section represent the distinct voltages at any node in the circuit. Respectively the flow variables represent the

1278

J.J. Granda / Mechatronics 12 (2002) 1271–1295

Fig. 7. Complete sensor model as entered in CAMP-G.

velocities and the currents involved. This approach considers the system as whole so that the state matrix involves all three sections of the sensor. The logic programmed in CAMP-G follows generalized matrix operations to obtain the desired transfer functions. Here, these are obtained in symbolic form using the A, B, C; D computergenerated symbolic system matrices. It generates the rate space form for the system equations, CAMP-G will generate the A, B, C, D matrices. The transfer function generalized matrix operations require the solution of the matrix form: 1

H ¼ CðsI  AÞ B þ D

ð10Þ

where H is the vector or matrix of transfer functions depending on how many inputs the system has. It is important to note that the approach proposed here can obtain several transfer functions simultaneously since for every row of the C matrix, a transfer function will be generated which relates that output to each input of the system so the approach is valid for single input, single output (SISO) or multi-input, multi-output systems (MIMO). In this paper, the computer-generated transfer function for the voltage across the capacitor crosses two different energy domains since the model is all together. The transfer function is obtained in one step in symbolic form. Any other transfer function for the efforts and flow output variables can be obtained. At this point the computer-generated model become so versatile, that all the linear control theory operations implemented in the control systems tool box can be used on the multidisciplinary model. For example, a bode plot can be generated using the computergenerated transfer function or the A, B, C, D matrices in order to do a Frequency

J.J. Granda / Mechatronics 12 (2002) 1271–1295

1279

Fig. 8. Computer-generated transfer function.

response analysis. Root locus, pole placement and other operations such as controllability and observability using the state space form are possible using the model produced by the approach presented on this paper. The result of the above matrix operations can be used directly in MATLAB or can also be generated in the form that SIMULINK uses to analyze such transfer function. The transfer function, which relates the output voltage to, the input force using CAMP-G is shown in Fig. 8. Reference [3] explains this in depth. There are other benefits of this approach. The computer-generated characteristic equation expressed in symbolic form as shown in Fig. 8 is written in terms of the physical parameters in symbolic form. This means that the coefficients of S are the result of the physical parameter combinations. This reveals the influence of each physical parameter on the coefficients and their dependencies. Such can be used for the study of sensitivities and most certainly to make changes to the model depending on design requirements. A computer simulation using this model was carried out using the symbolic transfer function and the state space matrices. The bode plot shown in Fig. 9 indicates the frequency range in which the sensor can measure the input precisely. The magnitude and phase angle are plotted showing a natural frequency of about 10,000 Hz as the reliable range of frequencies. This reveals that the

Fig. 9. Bode plot of sensor dynamics.

1280

J.J. Granda / Mechatronics 12 (2002) 1271–1295

Fig. 10. Step and impulse response.

sensor is accurate in frequencies below the natural frequency. Around this frequency, deviations in the magnitude as well as in the phase play a role in the accuracy of the measurements of the acceleration. The correct choice of the damping constants ensures that such deviations do not occur prior to 10,000 Hz. The lower limit is indicated by the pass band frequency. Using this model, different sensors of higher resonance can be designed, keeping in mid the influence of the piezoelectric crystal. Sensors capable of higher frequencies can be accomplished with small values of damping. Hoffman [6] studied this using a different method, which in comparison, verify the results of the one proposed here. Using the transfer function or the state space form, test functions were applied to study the transient response of the sensor. The step and impulse responses are shown in Fig. 10. The model produced from the bond graph and used in the computer logic for the automatic generation of the symbolic equations and matrices has been automatically transformed into the format that MATLAB uses to obtain these plots. In summary this demonstrates the theory and application of this proposed multienergy approach for the production of mechatronics models, as well as verified by the simulations in the frequency and the time domain respectively. 3. Electromechanical actuators Now attention is turned to another important aspect of mechatronics system, the electromecahnical actuators. Shown in Fig. 11 is an electromechanical model of a rotating actuator where the rotational position is to be controlled. Such model can be used to study position actuators or dc motors. Here the proposed process is that the PLANT which is the physical system, be modeled with bond graphs and the

J.J. Granda / Mechatronics 12 (2002) 1271–1295

1281

Fig. 11. Rotational actuator.

control system be modeled with block diagrams compatible with the SIMULINK notation. The reason is that the physical system is build with power bonds which represent the power distribution amongst the individual elements, this principle is valid for all systems in different energy domains and consequently allows the building of multidisciplinary models. The control part however follows signal flows rather than power bonds and not necessarily is subject to the basic physical principles of energy and power but rather to signal logic. Therefore the control part is more appropriately modeled with block diagrams. The purpose here is to demonstrate that in these combinations where multidisciplinary systems with controls are to be analyzed, how will the proposed methodology work? How can these be combined? How can one produce automatically the desired transfer functions? This paper illustrates some answers to these questions. The bond graph model consists of a mechanical and an electrical part (Fig. 12). The current is common to the armature resistance Ra and to the inductance La . This is represented by the 1 junction to whom both of the elements are attached. The inductance is represented by an ‘‘I-element’’, the resistance by an ‘‘R-element’’. The electromagnetic action is represented by an electric motor represented by the ‘‘GYelement’’, the friction and the rotor inertia are modeled by the R- and I-elements attached to the right 1 junction of the GY (motor). This 1 junction represents the angular velocity produced when a voltage is applied. The ‘‘C-element’’ represents the shaft compliance and is attached to a ‘‘0-junction’’ which is used to represent the relative velocity between the output of the armature and the actual mechanical inertia and rotational friction of the load. The proposed method here is that the physical system, the ‘‘PLANT’’ be modeled by bond graphs a general form that supports the mixture of mechanical and electrical components in an integrated model. Generalized variables, p and q (for momentums and positions on the mechanical side and electrical momentum (flux linkage) and charge on the electrical side are used as state variables. The output variables used are e, f (efforts and flows), e for forces, voltages and f for angular velocities and currents. Following the fundamental principles stated in [9] the state space form is found. The A, B, C, D matrices are

1282

J.J. Granda / Mechatronics 12 (2002) 1271–1295

Fig. 12. Bond graph model of rotational actuator.

computer generated in symbolic form using CAMP-G. The states are the momentums P12 and P7 , the position Q9 , and the momentum P3 . The computer-generated state vectors are Inputs vector u ¼ ½SE1 State variables vector p q ¼ ½P12; P7; Q9; P3; (arranged in logical order) 3.1. A and B matrices (by rows) Each computer-generated differential equation for each state variable derivative and the corresponding rows of the A and the B matrix are shown below: dP12 ¼ Q9=C9  P12=I12 R11 Að1; :Þ ¼ ½1=I12 R11; 0; 1=C9; 0; Bð1; :Þ ¼ ½0; dP7 ¼ P3=I3=G4x5  P7=I7 R6  Q9=C9 Að2; :Þ ¼ ½0; 1=I7 R6; 1=C9; 1=I3=G4x5; Bð2; :Þ ¼ ½0; dQ9 ¼ P7=I7  P12=I12 Að3; :Þ ¼ ½1=I12; 1=I7; 0; 0; Bð3; :Þ ¼ ½0; dP3 ¼ SE1  P3=I3 R2  P7=I7=G4x5 Að4; :Þ ¼ ½0; 1=I7=G4x5; 0; 1=I3 R2; Bð4; :Þ ¼ ½1;

ð11Þ

J.J. Granda / Mechatronics 12 (2002) 1271–1295

1283

The mathematical equivalence of the corresponding state space form of such equations then is 9 2 3 38 9 2 1 8 1 0 p12 > > dp12 > > > >  R 0 0 > > > > 7 6 7 > > 6 I12 11 > > C > > > > 9 dt > 6 7 7> 6 > > > > > > 7 6 > > 7 > > 6 > 1 1 1 7> dp7 > > > > p7 > > 6 7 6 > > >  0  R6 = 6 0 7 = 6 < < 7> I7 C9 I3 G45 7 6 7 6 dt þ 6 7 SF1 ¼6 ð12Þ 7 > > 1 1 6 7 7> > > 6 > > dq9 > > > > > 7 6 7 6 0 0 > 6 I > > > q9 > > 607 7> I7 dt > > > > 12 > > > > 6 7 7> 6 > > > > > > 4 > 4 5 > > > 1 1 5> > ; > : dp3 > ; 0  R2 : 0  1 p dt I7 G45 I3 3 Considering the velocity f12 as the output, the C and the D computer-generated matrices are % f12 ¼ P12=I12 Cð1; :Þ ¼ ½1=I12; 0; 0; 0; Dð1; :Þ ¼ ½0;

ð13Þ

Their state space form is  ff12 g ¼

1 I12

0 0

9 8 p12 > > > > = < p7 þ ½0fSF1 g 0 > > > q9 > ; : p3

ð14Þ

It is obvious at this point that the proposed multidisciplinary approach transforms the complete mechatronics system into the standard state space form f_xg ¼ ½Afxg þ ½Bfug and fyg ¼ ½Cfxg þ ½Dfug. Whether one starts with a bond graph or any other approach, this is a common ground for any method which will generate equations in second order form which then can be transformed to first order form. In this case the state space form is obtained directly from the initial bond graph model. It is obtained in terms of generalized variables, which intentionally are referenced by the bond numbers in order to make the approach general, but will represent variables across the disciplines. It is worth noticing that the proposed methodology is valid also for non-linear systems in which the state space form could be used to time varying coefficients, but in other non-linear applications. The computer-generated transformation of the bond graph into the Cauchy form of the equations will allow the simulation of a wide variety of non-linear multidisciplinary systems. This will be discussed in a subsequent section of this paper. What has been achieved with the proposed procedure is to generate the ‘‘PLANT’’ in state space form starting with a bond graph. CAMP-G generated the y vector of efforts and flows for the entire model and any desired transfer function in symbolic form. Since the .m files contain the A, B, C, D initialized matrices then the state space block in SIMULINK contains now the bond graph model as the ‘‘PLANT’’. This action in fact transfers the system’s bond graph to SIMULINK. The system entered in SIMULINK is shown in Fig. 13.

1284

J.J. Granda / Mechatronics 12 (2002) 1271–1295

Fig. 13. Bond graph state space model in SIMULINK.

Fig. 14. Controller output and input signal.

Considering an actuator designed to follow a pulse input representing a position control one can simulate the control system and the plant together. Fig. 14 shows the results of the model obtained using the bond graph model as a state space block in SIMULINK. Fig. 14 demonstrates the position tracking of the electromechanical actuator. The pulse wave is the input and shown is also the response of the system, which shows us how closely the controller follows the input signal. This simulation demonstrates that the proposed approach here can be used not only to analyze the physical system characteristics, but can be used to design control systems. It is obvious that at this point all the tools developed for this purpose in MATLAB and SIMULINK are available to design such systems. The approach in fact is a bridge between new bond graph modelling technique and the classical block diagram approach to represent control systems.

J.J. Granda / Mechatronics 12 (2002) 1271–1295

1285

4. Complex non-linear multidisciplinary models 4.1. Thermo-fluid, thermo-mechanical systems Complex non-linear multidisciplinary models play and important role in mechatronics. These can be analyzed with the proposed procedure as presented for the Sensors and Actuators with the appropriate variations that allow the representation of non-linearities. Models that involve thermo-mechanical processes and those that involve hydro-mechanical and hydro-electrical components are good examples. In order to illustrate this, this paper proposes a complex model of an internal combustion engine. The thermo-mechanical, hydro-mechanical models will be analyzed first and then a proposed thermo-fluid, thermo-mechanical-electronic model is proposed. Fig. 15 shows the combustion chamber with the mechanical components of the crankshaft and the connecting rod. The combustion chamber involves the thermodynamic process that transforms heat into mechanical power. It also contains the intake valves and the thermo-fluid process involved. There are modeling considerations necessary to obtain a closed set of differential equations. For example, the motion of the piston is geometrically constrained and dependent on that of the connecting rod and crank. This immediately suggests that if the linkages are considered to be rigid bodies in plane motion, the inertia effects are dependent on one another thus producing derivative causal relations. The masses of the piston, crankshaft, and connecting rod are statically dependent on each other’s

Fig. 15. Single cylinder internal combustion engine.

1286

J.J. Granda / Mechatronics 12 (2002) 1271–1295

position once they work together as a system. This fact does not allow the formulation of explicit differential equations of motion, a modeling problem that bond graphs detect right away, a problem that may cause computational problems when trying to solve an implicit set of differential equations. In order to remedy this situation, energy storing (C-elements) and dissipating components (R-elements) are considered. A field reduction is performed by lumping the inertial members into independent configurations. This allows to obtain an independent closed form set of differential equations. 4.1.1. Non-linear kinematic relationships The kinematics transformations between angular velocities and linear velocities as well as between torques and forces are described by means of non-linear kinematics relations. The symbol TF has been used to represent these relations. A constrained kinematic transformation occurs between the piston velocity and the angular velocity of the crankshaft. This relation controls also the transformation between the force applied on the piston and the generated torque on the crankshaft. Non-linear modulated transformers model these relations [1]. These elements and their relations are shown in Fig. 16. A transformer represents the interactions between lumped translation piston inertia and the rotating inertia of the flywheel. This ‘‘field reduction’’ requires finding the relations between the angular flywheel velocity and the piston velocity. Moreover, we must consider a field reduction of the inertial components. An external kinetic energy source is necessary to tell the simulation to keep the crank-slider in motion while the piston is near the top and bottom end of the stroke in order to establish a net forward momentum. A velocity source (SF28), based on the average speed of operation, was added at a velocity summation point (0-junction) as shown in Fig. 15. The modified bond graph model using field reduction techniques is shown in Fig. 15. Note that the reciprocating mass of the connecting rod has been included in the

Fig. 16. Piston-crank geometrical constraints and non-linear transformations.

J.J. Granda / Mechatronics 12 (2002) 1271–1295

1287

piston mass and that the rotating end of the rod and crankshaft have been included in the inertia of the flywheel. A harmonic balancer and torsional compliance has been added to the crankshaft. The harmonic balancer acts to dampen torsional vibrations of the crankshaft, which results in smoother engine operation while lowering fatigue loading on the crankshaft. Taylor [15] and Fergussen [16] discuss this topic further. Granda and Sime [13] explore the fundamentals of modeling the internal combustion engine cycle. Engja [17] explores diesel engines and bond graphs. Friction of the piston against cylinder walls and in the journal bearings has also been included in the final model. This model simulates the actual running condition of an engine. It produces a explicit set of independent differential equations. 4.1.2. Intake, exhaust ports and combustion chamber The bond graph model of the thermodynamic and thermo-fluid system is developed by first modeling the thermal and hydraulic phenomena [19]. They are integrated with the kinematics phenomena in the completed model. The combustion chamber is represented by a non-linear compliant element (a multiport C-element) that is represented by a pseudo-bond graph. . . (Fig. 15). The flow of the air is shown as a solid bond line (power bond) whose state variables are pressure and volume flow rate. The thermal energy contained in the air/fuel mixture is represented in the pseudo-bond graph (dotted lines) because the product of the state variables, temperature and heat flow rate, is not power and thus can not be represented as a true power bond. As a consequence the process of combustion has been separately modeled and computed by an iterative technique, which yields cylinder temperatures, pressures, and absolute energy content versus cylinder volume. This was done using the techniques outlined by Benson and Whitehouse [18]. These values have been computed by a separate program designed to analyze the relationships within the thermodynamic combustion process. These values then were made available to the simulation program in the form of tables in the input file. The combustion chamber is represented by an energy storing compliant element, a multiport C-element because the pressure increases at a rate proportional to the change in volume. The change in volume is the generalized displacement, which in turn generates an effort in a manner similar to Hooke’s law for a spring. Moreover, the pressure in a gas filled cylinder is also dependent on temperature hence requiring a third port on the multiport C-element (Fig. 15) on which temperature acts as an effort and the change in heat flow with respect to time is the generalized flow. The governing constitutive relations for the three port C-element during the intake and exhaust part of the cycle are based on ideal gas laws: Pressure ¼ ðR=CvÞ ðenergy=volumeÞ Temperature ¼ ðl=CvÞ ðenergy=massÞ

ð15Þ

The two-port resistive elements shown at the top of the Fig. 15 represent the intake and exhaust valves. The opening and closing of the valves can be modeled as a nonlinear resistor dependent on the motion of the crankshaft. A subprogram or function models the valve as an infinite resistance (zero flow) when closed and an ideal orifice

1288

J.J. Granda / Mechatronics 12 (2002) 1271–1295

governed by a compressible isentropic flow equation when open. The mass flow is given by the following equation [21]: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   A00 Cd Pu 2K ð16Þ DMDT ¼ pffiffiffiffiffi ðPRð2=KÞ  PRðKþ1Þ=K Þ RðK  1Þ Tu where A: valve area, Cd : discharge coefficient, Pu : upstream pressure, Tu : upstream temperature, PR: pressure ratio, K: universal gas constant. In addition, a control statement that checks for supersonic flow will fix the pressure ratio (PR ¼ PCR) between the port inlet and outlet when the state of choked flow is reached [15]. The critical pressure ratio is given by  ðK=ðKþ1ÞÞ 2 PR ¼ ð17Þ K þ1 The effects of valve leakage have been introduced by setting a very small constant mass flow rate during all parts of the cycle. 4.1.3. Thermo-fluid and kinematics systems The connection between the thermo-fluid and the kinematics subsystems is easily accomplished by the use of a transformer element, which converts the cylinder pressures into forces on the piston. This is modeled by a transformer element next to the multiport C-element. The transformer modulus, (piston head area), follows the mathematical relation between force and pressure. Since it is power conserving, it must do the same for the velocity and the volume flow rate as expressed in the relation: Force ¼ pressure  area Velocity  area ¼ volume flow rate

ð18Þ

Note that the mechanical motion of the piston changes the volume of the cylinder, which is one of the hydraulic state variables of the thermo-fluid system. In order to complete the link between these subsystems a source flow (SF27) is added which mathematically relates the rate of change of the cylinder volume to the crankshaft angular velocity. The source provides a volume flow rate due to the peculiar motion of a crank-slider driven piston. The flow is directed at the zero junctions on the pseudo-bond graph thus completing energy conservation [14]. All these non-linearities were implemented into the computer-generated model [20]. 4.2. Simulation results The graphs presented below demonstrate the operation of the engine and the dynamic analysis. The graph shown in Fig. 17 is the pressure vs. volume (PV) diagram for this internal combustion cycle. It is exactly the Otto cycle, as we would expect for a four-cycle engine. This result verifies our model as it can be compared to those shown [15], [18]. These can be obtained using the CAMP-G generated models and a simulator such as MATLAB or ACSL [12] to which it is capable of

J.J. Granda / Mechatronics 12 (2002) 1271–1295

1289

Fig. 17. PV diagram or fuel–air mixture.

Fig. 18. Cylinder pressure response in three cycles.

automatically delivering the models. The output of the simulation corresponds to experimental results and that of other methods of calculation. The cylinder pressure versus time output Fig. 18 is just like that observed on a digital scope display of a pressure transducer output of an actual running engine. This is usually the

1290

J.J. Granda / Mechatronics 12 (2002) 1271–1295

Fig. 19. Cylinder gas temperature.

experimental test result when observing the compression of each cylinder to determine the compression ratio and status of an existing engine. This model allows the evaluation of the temperature of the gases in the combustion chamber. The temperature has effects on the intake and exhaust flows. Fig. 19 shows the temperature change during the four strokes. The temperature is at its maximum at ignition during the power stroke while the volume is at its minimum. The temperature decreases rapidly during the exhaust stroke. The net force acting on the piston pin is shown in Fig. 20. It shows the influence of the forces generated by the kinematic links and the crankshaft on the pin. This graph is important in the design process, as these dynamic forces are difficult to measure.

Fig. 20. Net force on piston pin during three cycles.

J.J. Granda / Mechatronics 12 (2002) 1271–1295

1291

4.3. Eight-cylinder engine The single cylinder model can be extended to represent a complete eight-cylinder engine. This implies the extension of several single cylinders synchronized to simulate the continuous operation of the engine. It also means that the rest of the crankshaft model has to be extended to represent all the bearings, added inertia and stiffness effects. Using the graphics capabilities of CAMP-G an eight-cylinder model is produced as shown in Fig. 21. Granda and Channel [10] elaborate more on the modeling aspect of a complete engine. Since the original idea is to start with a bond graph and then transform the model into suitable computer models, it follows that other simulation programs tailored to non-linear systems can be used as the calculation engine in the simulation phase. In this case MATLAB will allow the programming of the non-linearities as separate functions, SIMULINK has the so called S-function capability to do this, ACSL uses macros. Other simulators such as, EASY5 and others that allow open programming which can also be used. The point is that the logic used in the proposed approach for a single one-cylinder engine is a complex thermomechanical, thermo-fluid system, which is flexible enough, be extended to a complete

Fig. 21. Eight-cylinder engine model.

1292

J.J. Granda / Mechatronics 12 (2002) 1271–1295

eight-cylinder engine. Such model is a complex MM-linear model, which using the methods presented here can be analyzed together with its electronic control system. In the above bond graph, the center structure represents the crankshaft with all its lobes and bearings. The symmetric structures on both sides of the crankshaft represent the cylinders. The computer model is generated in the same way as the more simple systems presented at the beginning of this paper. The non-linearities are entered in the input files of the simulator and a simulation of the complete engine performed. The author has further investigated these interdisciplinary subsystems. Mergen and Granda [11] presented the development of the valve subsystem using a combination of dynamic analysis and finite element analysis in order to study the valve springs. Fitsos and Granda [4] studied the development of alternative valve driving mechanisms, an alternative to the hydro mechanical systems to replace them with electromechanical rotational actuators. Granda and Channel [10] researched the modeling process for the eight-cylinder engine. It is clear that an integration of these models and the method proposed in this paper can produce very comprehensive models considering detailed internal parts of the engine. The method allows for a wide range of options. This paper concentrates on mechatronics systems as interdisciplinary models. In this light, it follows that if we were to consider the engine on its starting stage rather than the continuous operation, an interesting more complex system results if we model the engine with the starter solenoid and battery circuit. In that case we would have an electrical system combined with a thermo-fluid, thermo-mechanical system. In order to illustrate this principle, the starter solenoid and battery subsystem is added to that of Fig. 15. The resultant thermo-fluid-mechanical-electrical model is shown in Fig. 22.

Fig. 22. Internal combustion engine with starter.

J.J. Granda / Mechatronics 12 (2002) 1271–1295

1293

This model also is a close explicit model for which the method and process proposed in this paper is valid. This model joins all the challenges presented here. It is a multidisciplinary mechatronics system, it has several non-linear subsystems in the kinematic relations, as well as the thermo-fluid ones, and it requires external data for the representation of the combustion process, all in one. However the method proposed here produces a multidisciplinary model, which uses a generalized notation that allows the simulation and the integration of all subsystems as demonstrated from the basic models to the more complex ones throughout this paper.

5. Conclusions The paper has presented a wide view of what the role of bond graph modeling is in the analysis of mechatronics systems. It concentrates on three major components of mechatronics systems, sensors, actuators and complex physical (PLANT) models and clearly demonstrates the ability of bond graphs to be of assistance in all these areas. The ability not only to model physical systems in several energy domains or a combination of them has been demonstrated. The ability to translate such methods into logic algorithms that aid the engineer into the creation of very complex multidisciplinary models should contribute to bring theoretical considerations into applied methods useful in analysis and design, One step to which this research contributes is to answer the question on how to generate the physical systems models quickly and reliably. The models processed by the algorithms that transform individual elements into complex systems, check for implicit equations, algebraic loops and other mathematical problems that may arise simply by the way elements are connected. The ability to check for the correct causal relations and differential equations across different fields most certainly aids not only to make better models, but mathematically explicit correct models. This leads the engineer into considering the precise details of the physics involved and of course avoiding time consuming errors which otherwise will be difficult to detect. The methods proposed here also contribute to better communication between those who have used bond graphs in the past with those who are experts in block diagrams and use popular simulation programs such as MATLAB and SIMULINK. The software algorithms contained within the method presented here build that bridge and allow people who are trained on bond graphs to communicate well with those trained in the block diagram approach. This means they can work with a common language because the representation of systems in the MATLAB .m files for the state space form or the Cauchy form of the equations and system matrices can be well understood by both of them. Another advantage of this idea is the use of the same tools by people of different backgrounds and across the disciplines, fulfilling an objective of a mechatronics approach to design. This most certainly is one of the solutions sought by the experts in mechatronics, so that people from different backgrounds can work together with the understanding of the same models.

1294

J.J. Granda / Mechatronics 12 (2002) 1271–1295

References [1] Granda J. Modeling Methods for nonlinear discontinuities using computer generated models, MATLAB M-Functions and SIMULINK S-Functions. CIFA 2000 Conference Internationale Francophone de Automatique. Lille, France, July 2000. [2] Karnopp D, Margolis D, Rosenberg R. Systems dynamics: modeling and simulation of mechatronics systems. John Wiley & Sons; 2000. [3] Granda J. Computer generated transfer functions CAMP-G: interface to MATLAB and SIMULINK. Proceedings of the International Conference on bond graph modeling and simulation ICBGM’99, San Francisco, CA, January 1999. [4] Fitsos P, Granda J. Bond graph modeling of engine valve and control system using electromechanical rotary actuators. Proceedings of the International Conference on bond graph modeling and simulation ICBGM’99, San Francisco, CA, January 1999. p. 353. [6] Hoffman J. Matlab und Simulink Beispielorienterte Einfuhrung in die Simulation dynamischer Systeme. Addison Wesley; 1998. [7] Ashley S. Getting a hold on mechatronics. Mech Eng May, 1997. [8] Cellier FE. World Wide Web––the global library: a compendium of knowledge about bond graph research. Proceedings of the Third International Conference on Bond Graph Modelling and Simulation ICBGM’97, Phoenix, AZ, 1997. [9] Granda J, Reus J. New developments in bond graph modeling software tools: the computer aided modeling program CAMP-G and MATLAB. The IEEE International Conference on Systems, man, and Cybernetics, Orlando, FL, October 1997. [10] Granda J, Channel G. V-8 internal combustion engine bond graph model, a detailed modeling procedure. Proceedings of the International Conference on Bond Graph Modeling and Simulation ICBGM’97, Phoenix, AZ, January 1997. p. 233. [11] Mergen HJ, Granda J. Bond graph models and finite element models for engine valve spring transient response. Proceedings of the International Conference on Bond Graph Modeling and Simulation ICBGM’97, Phoenix, AZ, January 1997. p. 129. [12] Mitchel D, Gauthier J. Associates. Inc., Advanced Continuous Simulation Language (Acsl). User Guide/Reference Manual, 1997. [13] Granda JJ, Sime K. Computer aided modeling of the four stroke internal combustion engine. Proceedings of the 30th Heat Transfer and Fluid Mechanics Institute. California State University, Sacramento, June 1987. [14] Karnopp D. Modeling and simulation of adaptive vehicle suspension with pseudo bond graphs. Camp and Acsl. lMACS conference, Oslo, July 1986. [15] Taylor C. The internal combustion engine in theory and in practice. MIT Press; 1986. [16] Fergussen C. Internal combustion engines: applied thermosciences. New York: Wiley and Sons; 1986. [17] Engja H, Strand K. Modeling for transient performance simulation of diesel engines using bond graphs. Tokyo: ISME; 1983. [18] Benson R, Whitehouse. Internal combustion engines. New York: Pergamon N. Press. Inc.; 1979. [19] Karnopp D. State variables and pseudo bond graphs for compressible thermo-fluid systems. Trans ASME J Dyn Syst Meas Control 1979;101(3):201–4. [20] Karnopp D. Pseudo bond graphs for thermal energy transport. Trans ASME J Dyn Syst Meas Control 1978;100(3):165–9. [21] Domhold TG. Engineering experimentation. McGraw Hill, Inc.; 1966.

Dr. Jose J. Granda is a full Professor in the Mechanical Engineering Department at California State University in Sacramento California, U.S.A. He obtained his master degree at the University of California, Berkeley and his Ph.D. degree at the University of California, Davis. He is the Chairman of the Technical Activity Committee on Bond Graph Modelling of the Society for Computer Simulation International and a Registered Professional Engineer in the State of California. He has served also as General Chair or Program Chair of the International Conferences on Bond Graph Modelling and Simulation, ICBGM’93, ICBGM’95, ICBGM’97, ICBGM’99 and ICBGM’2001. During the 1999–2000 academic year he has been a visiting professor on sabbatical at the Bonn-Rhein-Sieg University of Applied

J.J. Granda / Mechatronics 12 (2002) 1271–1295

1295

Sciences in Sankt Augustin, Germany. Professor Granda has presented his research to universities and industries around the world since he is able to communicate in English, German, French and Spanish. He has lectured a the Ecole Centrale de Lille, France; (EPFL) Ecole Politecnic Federal de Lausanne, Switzerland, Bonn-Rhein-Sieg Teschnische Fachhoschule, San Augustin Germany, Technische Fachhochschule Trier, Germany, University of Louvan, Belgium; University of Guanajuato, Mexico, University of Zaragoza, Spain; National University, Colombia; Polytechnic University, Ecuador, Daimler-Chrysler, Germany, IBM Corporation. He has served as a consultant to industries such as LockheedMartin, ADAM OPEL AG, Germany, Lawrence Livermore Laboratories, Dataproducts Corp, and the FMC Corporation in the USA. Professor Granda has been a member of scientific and professional societies including the Society for Computer Simulation, ASME, ASM, Sigma XI the Scientific Research Society and the Phi Beta Delta Honor Society of International Scholars. In 2002, Professor Granda was awarded a NASA Faculty Fellowship. He worked at NASA Langley Research Center, Hampton, Virginia on the Dynamics and Controls of the International Space Station.