MSU LIBRARIES “ “i RETURNING MATERIALS: PIace in book drop to remove this checkout from your record. FINES will be charged if book is returned after the date stamped beIow. SINULAIION 0F NONLINEAR IULIIPORI ENGINEERING SYSTEIS By Guy E. Allen A THESIS Submitted to Michigan State University in pnrtiul fulfillment of the requirements for the degree of NASTER OF SCIENCE Department of lechunicel Engineering 1984 ABSTRACT SIMULATION OF NONLINEAR IULITPORI ENGINEERING SISTERS By Guy E. Allen A causal bond graph together uith its accompanying constitutive equations contains sufficient information to define the state equations of a given system. In the general nonlinear case, however. the problem of finding explicit state equations is often intractable. A.method has been developed to organize the system equations in a form suitable for numerical solution of the time response for a wide class of problems. This method has been incorporated in the bond graph modeling program ENFORI‘G and greatly reduces the time and effort needed to simulate many nonlinear systems. ACKNOWLEDGMENTS I would like to thank Dr. Ron Rosenberg for making this work conceivable; the faculty and staff of both the lechanical Engineering Department and the Case Center for Computer-aided Design for making it possible; and all my family and freinds for making it worthwhile. ii TIBLE OF CONTENTS LIST OF FIGURES 4.1 4.3 CHAPTER 5 - APPENDIX 1 APPENDIX 2 APPENDIX.3 APPENDIX 4 REFERENCES Introduction Analysis of Nonlinear Systems Fields in Bond Graph lodels Reduction of the Junction Structure Reduction of the Dissipation Field Reduction of the Storage Field Computational Considerations and Implementation Some Nonlinear Examples Coulomb Oscillator Rising-rate Springs Sprung Pendulum Conclusion An Example Run of the Nonlinear Bongraph Modeling Program ENPORTK.0 Calling Tree for ENFORTG.0 Subroutine Descriptions for ENFORTB.0 Source Code Listings for ENP0816.0 iii iv 13 14 15 22 22 24 29 32 34 43 45 69 LIST OF FIGURES Figure 2.1 Iultiport Field Structure Figure 2.2 Constitutive and Connective Relationships for lultiport Structures Figure 2.3 Junction Stucture Relationships Figure 3.1 Bond Graph Figure 3.2 Bond Graph After Augmentation Figure 3.3 Bond Graph Segment with Forced Causality Figure 4.1 Iechanical Oscillator Figure 4.2 Bond Graph of a Mechanical Oscillator Figure 4.3 Coulmb Damping Function Figure 4.4 Time Response of an Oscillator with Cbulomb Dumping Figure 4.5 Eighth Order Nechanical Oscillator Figure 4.6 Rising Rate Spring Deflection Curve Figure 4.7 Bond Graph of Oscillator with Linear and Nonlinear Springs Figure 4.8 Cbmparison of Linear and Nonlinear Response Figure 4.9 Comparison with Different Initial Conditions Figure 4.10 Pendulum on Spring Figure 4.11 Bond Graph for Pendulum on Spring Figure 4.12 A Time Response of a Pendulum on a Spring iv CHAPTER 1 Introduction Inherent in the design of a physical system or device is the prediction of its behavior. That is. given infonmation about the nature of the system it should be possible to model it and predict its function without constructing a prototype. This process of modeling. whether intuitive. physical or mathematical. is essential to the efficient development of a new design. Intuitive models have been used traditionally by craftsmen. builders. and the like and are still used in situations where the cost of the product is low and volume does not warrant additional expense. or where the analytical grasp of the underlying phenomena is poor. Physical modeling has often been useful where the intuitive understanding is insufficient and the mathematical model is intractable. The use of physical modeling is often limited by the time and expense involved. Advances in the physical sciences have increased the range of problems that can be approached computationally. Ihile the amount of calculation necessary to predict the behavior of even a relatively simple system can be prodigious. the ongoing reduction of the cost of high-performance digital computers has made computer-based modeling a useful and highly accessible design tool. Typical simulation categories are finite element. discrete event. and continuous time modeling. This work addresses itself to the last of these three. ie; systems which can be described by sets of ordinary differential and algebraic equations. There are a number of commonly used techniques which allow the user to describe the system to the simulation software package. One method is for the user to code the differential equations directly into FORTRAN or some other high level language and then have this code compiled and loaded with a main program. An example of this is the program called DIFFEQ. which was written at the A. H. Case Center for Computer-aided Design at Iichigan State University [1]. DIFFEQ consists of a command processor. integration routines. a postprocessor. and a user supplied subroutine which defines the model. This method executes quickly after compilation and loading but places the burden of finding the state equations and successfully translating them into FORTRAN on the user. which can be a considerable task for large systems. especially for the engineer who is not practiced in FORTRAN coding. Other approaches require the user to reduce the system to some analytical form. Programs such as CSNP. CSL and ACSL will accept relatively free-form input of either the state equations or the blockrdiagram-style transfer functions. and allow the simulation of virtually any system that can be described in this fashion. These are powerful packages and are useful in many applications [2]. Some modeling packages are tailored for particular kinds of physical systems and allow the description of the model in some physically intuitive form. SPICE. for example. is a circuit analysis program that allows the description of a physical device in terms that are natural to a circuit designer [3]. This makes the program very easy to use. but generally limits its scope to electrical problems. Other specialised softvare packages address other physical regimes. and more software of this type is becoming available. A useful technique for modeling more varied systems is the bond graph method. This allows the representation and manipulation of multiple-energy-domain systems in a symbolic form and provides a straightforward procedure for producing the system equations of the model [4]. Some of the software packages which accept bond graphs directly include ENFORT. TUTSII [5] and CARP [6]. oENPORT is a self-contained package written in FORTRAN. but current versions are limited to linear systems. CARP is a bond graph preprocessor that deals with nonlinear problems but cannot handle implicitly defined resistive effects or dependent storage effects and requires a simulation package such as CSL for the actual integration and post-processing. The work presented here deve10ps a general approach to computation of a numerical value for the time derivative of the state vector at a particular point in state space when given a causal bond graph. Chapter 2 contains the problem description and analysis. Chapter 3 discusses some considerations relevant to implementation in software and Chapter 4 has some simple examples to demonstrate the procedure. Chapter 5 presents conclusions and suggestions for further work. Chapter 2 Analysis of Nonlinear Systems 2.1 Fields in Bond Graph Iodels Nany physical systems which can be modeled by bond graphs can be represented as in figure 2.1 [7]. The system inputs are defined by the collection Of 3. and 8: elements and are referred to as the source field. Dissipative effects from the R elements of the bond graph are modeled in the loss field. Dynamic effects are the result of the inclusion of C and I elements in the model and are represented in the independent energy-storage field. If the problem has derivative causality it is also necessary to include the effects of the dependent energy variables. These are represented in the dependent storage field. Each of these fields has constitutive equations associated with it. and together with the connective structure represented by the junction structure fields. these define and interactions of the different components of the system. Specifically. then. we have: 0 - flu) (2.1) where the input vector U is a function of the time. 5 Source Field Se» Si _LT Independent C --l Paynter Storage I l— Junction Structure Dependent C l— l—‘cy Structure Storage l ——l (O, l) ,__ 1% Dissipation Field “4 TF l-— Modulated Junction Figure 2.1 Nultiport Field Structure The constitutive relations for the energy storage elements are zi ' ‘1‘11'xd) (2.2.) zd ' ‘d‘xi'xd’ (2.2s) These define the coenergy vectors. Z1 and zd, .3 function, of the current state. The i and d denote the independent and the dependent parts of the storage field. The output vector of the loss field (Do) is 3 function of the input vector to the 1088 11014 (Di). This relationship may be described by no . ‘1‘Di’ The connective relationships (2.3) represented by the 0. 1. TF. and GT elements can be dealt with as shown in figure 2.2. Across the modulated junction structure we have: To ‘ ‘t‘xi'fl‘Ti The invariant connective junction field: V0. :3 Ptvi where r - zi xd Vi 3 ‘ Do i U To . a and - 1 xi V0 - 2d (2.4) structure is represented in the Paynter (2.5) (2.6) (2.7) moumuoouum «wounds: new name-3330a 2530.330 mm: 0:332:00 «.N camar— a a Allin .lfl [News a a .N. l... A. lwg [II .x 2.2 Reduction of the junction structure. Every bond is classified as a unique type. implying a defined role for each effort and flow variable. From figure 2.2 it can be observed that bonds incident to (3., 8:. R. I. C. TF. G!) should also be incident to a 0 or 1 junction. Such bonds are classified as external. Their power Vttilb1°8 39904! in thO Vi and V0 vectors as shown in (2.6) and (2.7). Bonds incident to only 0 and 1 junctions are classed as internal. Their variables are collected in a single vector. vint' gnd include the internal effort variables and the internal flow variables. The Paynter junction structure is the one field which is invariant for any bond graph because it is composed only of 0 and 1 elements. Examining this field we may write vo ' Joi I vi * Jon ‘ vint (2.8a) vint ' J'ni . vin + Jnn . vint - (2.8b) Since the coefficients of J are all constants we may find a solution for vint from vint - ((1 - In)“ ”um“ (2.9.) 01' vint " 3.31“. (2.91:) It is assumed the the required inverse exists. Otherwise the junction structure is degenerate and no reduction of equations is possible. Using (2.9b) in (2.8a) we can find an expression for the reduced junction structure vo ‘ ”oi + Jon‘Jii"V1 (2.10) or v0 3 Joisvi (2.11) Figure 2.3 shows the role of the modulated junction structure. The equations are derived from (2.5) and figure 2.3 as follows: ‘p P 'I F q - - X Iii Is: Ii: 314 21 1;; zd . J” I,“ I“ J“ id + Is: [To] ' (2.12) ”1 131 Iss J'ss Isa Do 1,, . - . and _ _ '2: [13]" [In In In 3:4] ii + [In] [To] (2.13) Do P. 10 M Paynter _T_.> Modulated Junction Junction v' Structure T Structure 0 <9— (0, 1)_ J‘— (TF, or) Figure 2.3 Junction Structure Relationships The output vector V is not of analytical significance. and has been drOpped from the equations above. More succinctly: V3 ' P11.v; + Psi‘To (2.14) Ti - puev'i + Pasta (2.15) where xii Pi z; ' 0 Do U 11 and the matrices Pxi. P31. P13 and P3, are the aggregates of the matrices in (2.12) and (2.13). Combining the modulated junction structure field equation. (2.4). with (2.15) we can find T1 = ((I - 9,,tet)“ 'P..)svi (2.13) 01' 1'1 . pvtv; (2.19) provided the indicated inverse exists. If the modulated junction structure varies with the time or the state then the evaluation of ‘t and the reduction for the junction structure must be repeated during the calculation of each step of the time response of the system. Combining (2.19). (2.4). (2.14) and (2.15) we can arrive at the expression for the complete reduced junction structure: V3 ' (P11 + Psi ‘t P')V; (2.20) Shifting notation slightly by expanding v; and V3 into subvectors and partitioning the coefficient matrix in (2.20). we can now consider the other subfields of the system 12 3: F1811 Isis J'31s 1814- "z: zd ' Jssi Isa: 18,, 1334 id (2.21) Pg ‘18“ IS“ 18,, IS".‘ Do 3’. where for the most general case each of the elements of IS is a function of the time and state. Ihen defining causality using the normal sequential procedure the input to the derivative storage field (2d) 1. purely . function of the independent 3t01480 field (21) and system inputs (U). Therefore 18,, and 18,, _are zero. 18,, is also more because of the skew-symmetry of any power conserving junction structure. This allows the simplification of (2.21) to Xi - Is“. 21 + Jsueidrr 18,900 + 18,,‘0 (2.22) z‘l - Is“. 21 + J'S,,‘U (2.23) D, - :3“. 21 + JS,,'D° + 13...!) (2.24) To evaluate the time response of the system it is necessary to describe the time derivative of the state vector in some computational form. For the general linear time-invariant case it is possible to solve for the time derivative explicitly: i.e.. Xi . Aex1 + B‘U (2.25) 13 where A and B are constant-coefficient matrices. The nonlinear problem does not so readily admit to an explicit reduction. The problem may allow the symbolic reduction to a form equivalent to (2.25) where the A and 3 matrices are functions of the state and time. but this is not currently considered a practical approach for general computer based solutions. While the development of the explicit state equations is not always practical. the equivalent information is available in the form of the set of constitutive connective relationships discussed previously. These allow the calculation of the time derivatives of the state vectors for any particular position in state space. 2.3 Reduction of the dissipation field. Combining (2.4) and (2.24) we can get an expression for the input to the dissipation field: Di ' JS,,' 21 + 18,,‘61(Di) + 18,,‘U (2.26) For a system with explicit causal erields 18,, .111 be zero. gnd equation (2.26) will be an explicit expression for the Di vector. For a system with causally implicit dissipation fields the expression will be an implicit algebraic relationship for Di° The 1gttgr sitnmtion can be dealt with in a number of ways and techniques suited to bond graphs are discussed by Graf [8]. Having obtained Di' by 'hgtever means. no cgn then be obtained directly from equation (2.4). 14 2.4 Reduction of the storage field. At this point enough information is available to formulate an expression for 31- Exmmining equation (2.22). however. it is apparent that the formulation for 21 requires an evaluation of Id. To facilitate this equations (2.2a). (2.2b) and (2.23) may be combined to yield: 6d(xi. xd) - 18,,‘ditxi. 1,) + 18,,‘U (2.27) Thi‘ ilP11°8 4 solution for 1d and hence for 21. The result may be summarized in the form of a system of coupled implicit algebraic-differential equations of the form it, - 2(21. id. 1),. n) (2.28) 0 - “1,. X... a) (2.29) Vith 2, given by (2.2a) and Do given by (2.4). This set of equations then may be integrated numerically to calculate the system time response [9]. From the point of view of nonlinear bond graph models the most difficult reduction problems arise from implicit nonlinear dissipation fields and from nonlinear dependent storage fields. The reduction plan presented in this chapter should provide a guide to the possible reduction available in a specific model. based on an analysis of field characteristics. CHAPTER 3 Computational Considerations and Implementation The procedure outlined in the preceding chapter is applicable to a very broad range of problems. It is capable of taking almost any system which can be described by a bond graph and computing the time derivative of .the state vector for any state and time. This procedure has been implemented with some limitations in a nonlinear version of ENPORT. To compute a junction structure relationship of the type expressed in (2.21) it is necessary to classify and organize a good deal of information about the bonds that join the various elements of the system. Ihen doing this for a linear system there need not be any distinction made between bonds incident to transformer-gyrator elements and bonds incident to 0 and 1 junction elements. 'hen dealing with the general nonlinear problem using the procedure given in Chapter 2 it is necessary to classify these bonds into a number of distinct categories: 1) Bonds which join a field element and the junction structure. 2) Bonds which join two invariant junction elements. (vint) 3) Bonds which join modulated junction elements (TF. GT) to a field element or an invariant junction element. (Ti' To) 15 16 4) Bonds which join transformer-gyrator elements to other modulated junction elements. (Tint) The nonlinear version. ENPORT6.0. is based largely on the existing linear modeling software and does not have access to all of the information described above. The implementation in ENPORT6.0 allows the description of a very general class of connective subfields. but the bond graph must be augmented by the addition of two-port O and 1 junctions to separate transformer-gyrator elements from _both field elements and other transformer-gyrator elements. For example. a system described by the graph in figure 3.1 would have to be converted for ENFORT6.0 as in figure 3.2. This augmentation causes no loss of generality. Se-eTF-Hl—HGY—ul C R Figure 3.1 Sample Bond Graph Se —-li-*lTF—’-ll/-*-lGY—‘li—-l| \ C R Figure 3.2 Sample Bond Graph After Augmentation 17 Equations (2.2a) and (2.2b) are capable of describing an extremely broad range of storage elements. including causally mixed coupled multiports. Only single-port field elements are within the scope of the present software. This is primarily due to the difficulty of describing the complex functional dependencies inherent to this kind of system in a simple interactive fashion. The current capability for describing the constitutive relationships in a fashion compatible with system causality is based on the work of T. T. Chou [10]. The user may describe a particular functional relationship with whatever causality is convenient. The software will then check to ensure that the form required in equations (2.2a) and (2.2b) may be found. For example. consider the section of bond graph shown in figure 3.3. The causality on the R element requires a constitutive relationship such that the flow on the connecting bond is a function of the effort on that bond. If the user elects to define the relationship in terms of flow as a function of effort the software will determine whether or not the required functional inverse exists. It will prompt the user to use a more suitable function should that inverse not be available. R Se—uO—“l L Figure 3.3 Bond Graph Segment With Forced Causality 18 Once given the connective relationships of the junction structure and the constitutive relationships of the field elements. the time derivatives of the state can be evaluated as in (2.26) through (2.29). There are two considerations which give rise to computational difficulty: 1) Evaluation of implicitly defined nonlinear resistive (loss) fields. 2) Integration of the implicit differential equations which are caused by nonlinear dependent storage effects. These traditionally have been difficult problems for general purpose software of this type and are clearly delineated in the sequence of steps leading to the calculation of a time response. The current version of the nonlinear package. ENFORT6.0. can deal only with an explicit loss field and pure integral causality. However. the initial design was made with the most general purpose model in mind. Places where these difficult problems arise are clearly marked in the software. and improvements should be implementable without undue difficulty. The implicit algebraic problem in 1) above can be approached using traditional methods. This might not present an unacceptably large computational burden for relatively well behaved systems because a 19 'good' initial guess for the solution for the algebraic problem can usually be found from the state at the previous time step. Unfortunately. this may not be true for systems with large or frequent discontinuities. Another method which has promise uses knowledge of the system itself by creating a new dynamic problem which has a steady-state solution equal to the unknown state of the implicit static problem [8]. The primary advantage of this is the control of the problems of non-unique solutions. convergence. and stability often found when trying to solve a general nonlinear set of implicit algebraic equations. The combination of dynamic augmentation and more traditional methods may prove to be an effective treatment of the problem. The implicit integration problem poses other difficulties. Equations (2.28) and (2.29) form a general set of coupled nonlinear differential-algebraic equations of the form A(y.t) i - G(y.t) (3.1) For the case of integral causality the A matrix above will be an identity matrix and (3.1) defines a simple set of explicit coupled differential equations. Computation of a time response poses no problem in this case. Derivative causality will produce an implicit form. however. Examining (2.21) a general relationship in the form of (3.1) may be found: 20 i 18,, 18,, 18,, 18,, 2, 2d 18,, o o 18,, 1d (3.2) The number of equations available when (3.2) is expanded is equal to the number of state variables. The total number of unknowns must also include the number of dependent storage variables. The additional constraints required come from (2.27). which is a set of implicit algebraic constraints equal in number to the dependent storage variables. Thking (3.2) and (2.27) together produces a set of implicit coupled algebraic-differential equations. The calculation of the time response of such a system generally requires that the derivatives be known for the initial condition. One possible approach begins with differentiating (2.27) with respect to time: as . ad . a18 as . 318 ad . —-¢x+—-1-x-——’-16+18—1x+__'14+18—1x+ axi ‘ axd d axi i 31ax1 i axd ‘ uaxd d aJs “a 18,,1 i. ”8“!) + Isufl it, + axi ax, 1 ax, ax, aa ' Grouping terms of Xi and id we get 21 ad 318 ad 618 3U __g _ si‘i _ 18,, 1 _ sen _ 18,,—- xd a axd axd axd axd axd 34, 318,, ad, a18 . an . " ——‘1 - 18,, - ———’-—U "' 18,,—- X, (3.4) If the necessary inverses can be found then Rd can be obtained from the above. Using id in (2.22) produces X, directly. This procedure to find id may restrict the choice for initial conditions to those which allow the required derivatives and inverses to exist. CHAPTER 4 Some Nonlinear Examples 4.1 Coulomb oscillator. The system shown in figure 4.1 is a simple mechanical oscillator. The corresponding bond graph is shown in figure 4.2. In this example the mass and spring are linear. which is representative of many mechanical devices. while the resistance due to friction forces will be the nonlinear function shown in figure 4.3. This is a model of Coulomb damping. or dry friction. \\\\\\\ Figure 4.1 Mechanical Oscillator R l CAN—fin Figure 4.2 Bond Graph of a Mechanical Oscillator 22 23 F orce Velocity Figure 4.3 Coulomb Damping Function The discontinuity of the friction force allows for non-unique equilibrium values. The resistive force exerted on the moving mass is constant in magnitude. but changes sign with the direction of the motion. The force exerted on the mass by the spring is a linear function of the displacement of the spring from its relaxed state. For some sufficiently large initial condition the system will oscillate until such time as the spring force is no longer large enough to overcome the friction force. Any position of the mass where the force from the spring is less than force available from friction effects is a possible equilibrium point. The procedure to describe this system to ENFORT6.0 and get a time response for a particular set of initial conditions may be found in Appendix 1. The graphic results for particular initial conditions are shown in figure 24 4.4. In this example the system is released with some initial displace-ant (1,) and zero velocity (x,). The final equilibrium point for this initial value problem is clearly non-zero. and can vary with system paramaters or initial conditions. 2.00 '" v SCALE \ elOE ) —\ one. ,\ W ~l.20 — .2.00 n 1 1 1 1 1 1 1 l 0.00 0.80 1.20 1.80 2.40 3.00 TIME SCALE 010E 1 LEGEND: X01-—-X02 mm Figure 4.4 Tnme Response of an Oscilator with Coulomb Dumping 4.2 Rising-rate springs. A slightly more complicated example. as shown in figure 4.5. allows the comparison of a linear and a nonlinear storage field. In this system the springs connected to the lower two masses are linear. while the upper springs are modeled as risingnrate. These springs have a force-deflection curve which is a cubic polynomial. as shown in figure 4.6. A system with 25 . Ia. \\\\\\\\\\\\ \\ Figure 4.5 Eighth Order Mechanical Oscillator F orce Displacement _ Figure 4.6 Rising-rate Spring Deflection Curve 26 these springs may exhibit amplitude-sensitive frequency response. as motion with large excursions may be driven more forcefully than motion near zero displacement. 1 1 XE A“: R CR (3 Figure 4.7 Bond Graph of Oscillator with Linear and Nonlinear Springs The bond graph in figure 4.7 generates the independent state vector 27 l- q ‘1. Pa pas xi ' qas qis Pia psa qss where Q, and q,, are the displacements on the rising-rate springs. and q,. and 0,, are the displacements on the linear springs. The p variables. represent the momentum of the masses. The subcripts refer to the bond numbers in figure 4.7. The response shown in figure 4.8 is for identical initial displacements for both the linear and nonlinear springs and light damping. Figure 4.9 shows the response for a system with initial conditions equal to twice that demonstrated in Figure 4.8. Comparing figures 4.8 and 4.9 marked differences in the nonlinear subsystems for each simulation. Note that as well as the differences in the detail of the response of the nonlinear system for the different initial conditions there is also a marked difference in the period of the motion. The linear portion of the system shows only a change in the magnitude of its motion. as expected. 28 1.50 T Y SCALE s10E 0 0.90 0.30 -0.30 -0.90 -1.50 1 1 1 1 1 1 1 1 1 0.00 0.30 0.50 0.90 1.20 1.50 TIME SCALE .10E 1 LEGEND: X01——-X04 mm XOS-u— X08-—— Figure 4.8 Comparison of Linear and Nonlinear Time Response 3.00 Y SCALE 0105 0 1.80 0.60 -0.60 -1.80 '3.00 ' l l l l l l 1 1 0.00 0.30 0.60 0.90 1.20 1.50 TINE SCALE O10E 1 LEGEND: X01— X04 --- X05 - X08 ..... Figure 4.9 Comparison with Different Initial Conditions 29 4.3 Sprung Pendulum. Figure 4.10 is a mechanics problem with a nonlinear connective field. The appropriate bond graph is shown in figure 4.11. The velocity relationships across the junction structure are r - cos(0)Vx + sin(G)V& b . _ESM v + MY r ‘ r y These correspond to the port equations - cos(0) f, 8 sin(0) f, , - cos(0) f, ,, - -sin(0) f,, as ' r/£,, “Q f f f f f Figure 4.10 Pendulum on Spring 30 Figure 4.11 Bond Graph for Pendulum on Spring The state variables of interest are the angle of the pendulum relative to the horizontal (1.) and the radius (1,). A typical response of this system is shown in figure 4.12. and a coupling phenomena between the rotatinal and axial motion can be observed. 2. 50 Y SCALE 110E 2. LEGEND: l 28 .06 .84 .62 .40 0.40 X02——-X04.mm Figure 4.12 31 L 1 I l l 1 l 0.80 1.20 1.60 TIME SCALE tIOE 1 Time Response of a Pendulum on a Spring 2. 00 Chapter 5 Conclusion This work allows the calculation of the time response of a nonlinear physical system represented by a bond graph. Starting with the information available from a causal bond graph the derivatives of the state vector are calculated at a given state and time. This allows the subsequent computation of a trajectory through state-space. The current software implementation can interpret any system containing one-port field elements (C.I.R) and two-port junction structure elements (TF.GT) which has integral causality and explicit resistive fields. The general analysis presented in Chapter 2 outlines a procedure without these restrictions and can deal with causally coupled multi-port field elements. derivative causality. and implicit resistive fields. The underlying design and implementation of the software were done with the most general case in mind. and further development of the program should be straightforward. In use the software is similar to ENPORT5.2. an interactive linear bond graph modeling package. This package is a powerful interactive modeling tool. and the addition of nonlinear capabilities has greatly broadened the scope of problems with which ENFORT can deal. Further work which is indicated includes: 32 33 1) Completion of the implicit resistance problem. 2) Completion of the derivative causality problem. Improvements to the software may include: 1) 2) 3) 4) Nonlinear multi-port elements Implementation of the implicit resistance and derivative causality problems. The ability to make a modulus a function of any variable in the system. Modulated field elements (R. I. C). APPENDI CIES Appendix 1 An Example Run of the Nonlinear Bond Graph Modeling Program ENPORT6.0 This is a copy of the interactive pregram run which generated the example listed in chapter 4. 'adda ya want now??? RUN77 ENFORT6.0.TEST O0.0000000000000000000000000000000000000 00.0....0.0.0.0..OOOOOOOOOOOOOOOOOOOOOOO ee es “ ENPORT‘6 “ ee ee 9‘ BOND GRAPH PROCESSOR “ 9‘ FOR.NON-LINEAR SYSTEMS ‘8 ee ee “ VERSION 0.0 ‘9 es ea ‘8 FOR INFORMATION CONTACT “ “ DR. R.C. ROSENBERG “ as ea OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO OOOOOOOOO0.00000000000000000000000000000 PRESS (RETURN) KEY TO CONTINUE... WHAT KIND OF TERIINAL ARE YOU USING? : 4006. 4010. 4014 OR OTHER STORAGE TUBE DEVICE 4114 4025 OR 4027 : 4105. 4107 OR 4109 SOMETHING ELSE NUOU> ENTER CHOICE(D):D INNT TO RELOAD A STORED PROBLEM? (Y):N 34 first 35 IANT TO EITER A NEW PRQLEM? (Y): PRQLEN TITLE OPTION 8 T: SEE OR SET THE PRWLEI TITLE (1 LINE) L: LOOK AT THE PRmLEM DESCRIPTION N: ENTER A NEW DESCRIPTION C: CHANGE AN EXISTING LINE H: HELP (=DEFAULT) X: EXIT FROM me PRQLEM TITLER DITER OPTION (H):X no LEAVING 1113 manual: 1mm «a PROCEED? (Y) : BQID GRAPH WORKSPACE MANAGER OPTIONS INPUT 1) NEW now am anrr me CURRENT BOND am LOOK AT ma saw) am 881‘ on mu ma POWER museums DISPLAY ma POWER mascnons mnsm 1111: cam 10 SYSTEM FOR mocessma DEFAULT sumo OPTION CHANGER mp axrr FROM wompaca manna (murmur) filgflg'fil‘tfllfl mm OPTION (I):I THE 'ORRSPACE I8 EMPTY. 'ANT INSTRUCTIONS? (N): MEEGRAHMDIHONS: >11 .C2.R3 >1123 ) e» cam IDDIFICATIONS comm «a '8 NANAGER OPTIONS (I.E.L.P.D.T.N.H.X. ):L NODE NODE NEH NAME BWDS 36 0) Haunt-4 HWNH HIT (RETURN) TO (DNTINUE... WS MANAGE OPTIONS (I.E.L.P.D.T.N.H.X.):T eee PREPARING TO EANSFER THE W8 GRAPH... POWE DIRECTIONS SHWLD BE SET BEFORE TRANSFER. IF YOU DID SO. CONTINUE. OTHERWISE. SET POWE DIRECTIONS. CWTINUE? (Y): see WORRSPACE GRAPH HAS BEEN HANSFERED TO 'THE SYSTEM. WS MANAGE OPTIONS (I.E.L.P.D.T.N.H.X. ):X PROCEED? (Y) : FUNCTION OPTION 8 LIST CURRET FUNCTIONS CHANGE SEECTED FUNCTIONS SET ALL FUNCTIONS UNIFORM FUNCTIONS (ZERO OR (NE) HEP EXIT FUNCTIONS SECTION (IDEFAULT) 5589““ ETE FUNCTION OPTION (L. C. S.U.H.X. (FULL)):L 1:E2=LINEAR (Q2) 2:F1'LINEAR (P1) 3:E3=LINEAR (F3) _..._w _ 38 HIT (RETURN) TO CONTINUE... BITE FUNCTION OPTION (L.C.S.U.H.X.):C WHAT FUNCTION DO YOU WANT TO CHANGE? 3 FUNCTION 3 NOW DEFINE AS E 3 I LINEAR (F 3) SET FUNCTION TYPE (LINEAR):COUL E 3 . C1 ‘ SIGN(F 3) VALUE OF C1 ? ( 1.0000 ): VEIFY ? (NO): eee gvs _ 0 see REVERSE INPUT AND OUTPUT 7 (N0): m FUNCTION ORION (LsCsSsUsneXe (FULL)):X PROCE-I (I): P‘. INPUT DEFINITIONS “‘ DO YOU WANT HELP? (N): SEVE OPTIONS : INITIAL CONDITIONS TIME comms : SEVE THE PROBLEM SAVE THE RESULTS IN A FILE HEP : EXIT THE SEVE (I'DEFAULT) “saint-3H ETER OPTION (X):I INITIAL CONDITIONS 8. SET TEE INITIAL CONDITIONS P: PRINT THE INITIAL CONDITIONS C CHANGE SEECTED INITIAL CONDITION 8 Z SET ALL INITIAL CONDITIONS TO ZERO OR (NE 39 R: RETURN (IIDEFAULT) ETER OPTION (R):S SET THE INITIAL CONDITIONS: X 1 ? ( 0.0000E-01) : X 2 ? ( 0.0000E-O1) :17 ETE INITIAL CONDITION OPTION (S.P.C.Z.R.(FULL>): INITIAL CONDITIONS SET THE INITIAL CONDITIONS PRINT THE INITIAL CONDITIONS (muse SELECT- INITIAL CONDITIONS SET ALL INITIAL CONDITIONS TO ZEO OR (NE RETURN (I'DEFAULT) EEQ‘” ETER OPTION (R): SEVE OPTIONS INITIAL CONDITIONS TIME (DNTROLS SEVE me PROBLE SAVE THE RESULTS IN A FILE HELP : EXIT THE SEVER ('DEFAULT) “539?? BITE OPTION (X) :T TIME CONTROL PARAMETES FOR SOLUTION: ETE INITIAL TIME ( 0.0000E-01) : ETER FINAL TIME ( 0.0000E-01) :30 ETER NUDE OF POINTS TO STORE ( 2) :51 SEVE OPTIONS INITIAL CONDITIONS TIME WNTROLS SEVE THE PRmLEM SAVE THE RESULTS IN A FILE : HEP : EXIT THE SEVER (-DEFAULT) Hlfl'ldmz-IH ETER OPTION (X):8 SELECT (MITPUT VARIABLES FOR LATER DISPLAY. OUTPUT VARIABLES SEECTION L: LIST THE CURRET VARIABLES 8: SET THE LIST 40 A: PUTALLX.UORDON'mELIST D: DROP ALL X. U OR D FROI THE LIST R: RETURN (-DEFAULT) mm OPTION (R): WRANR OPTIONS A: ADAMS-BASEFOUR’IB PRDI CTOR-CORRECTOR INTHERATION BACX'ARDS-DIFFERE‘ITIATION NE'TBOD OTHER EXOTIC ME'mODS RUME-KUTTA (THIRD-ORDER) RETURN (=DEFAUL'T) 959E m OPTION (Q):R 5.25758 CPU seconds used. mm OPTION OR HIT (RET): WRATOR OPTIONS ADAMS-BASEFOUR’IB PREDICTOR-CORRECTOR INTEGRATION BACXIARDS-DIFFERENTIATION METHOD OTHER EXOTIC METHODS RUME-KUTTA (THIRD-ORDER) RETURN (IDEFAULT) 9599? m OPTION (0): ~SCLVER OPTIONS INITIAL CONDITIONS : TIRE (DNTROLS SGNE 1m: panama SAVE IRE RESULTS IN A FILE : HELP EXIT IRE SGNER ('DEFAULT) rflFmF-ig-l mm OPTION (X) :X PROCEED? (Y) : POST-PROCESSOR OPTIONS DRAW A GRAPE PRINT A TABLE PRINT A PLOT HELP EXIT POST-PROCESSOR (DEFAULT) 55??? DITER OPTION (X) :P 41 EB CURREVT DISPLAY TITLE IS: DUIIY DISPLAY LABH. - CHANGE IE mum 100 um 10 CHANGE rr: (Nmr mass m NE! nun ON ONE LINE: comma DAMPER me me comm DISPLAY nun IS: comma mum mini; mum 200 um 10 CHANGE rr: (N): ADD DATE AND um (I): (230083 ma DISPLAY VARIABLES-— ml: cum msw LIST: «t m. ANY ADDITIONS? (N):! NEXT VARIABLE? (BLANK To 0011):: 1 NEXT VARIABLE? (Bum: m (1mm: 2 mm VARIABLE? (BLANK 'm QUIT): ma cum DISPLAY LIST: X 1 X 2 ANY DELETIONS? (N): ANY ADDITIONS? (N): SET TIIE LINITS FOR DISPLAY: mm INITIAL TIIE ( 0.0000E-01) : mm FINAL TIRE ( 3.0000E+01) : SET TINE INGENENT FOR DISPLAY: mm DISPLAY INCRENENT ( 6.0000E-01) CONTINUE? (Y): 42 08/09/84 20:22:42 COULOIB DAMPER EXANPLE PLOT-‘-X 1 PLOT-O-X 2 0.00E-01 6.00300. 1.20301 1.8080-01 2.4OE+01 3.00EI-01 1.35301 0 9 . . . . I ‘ 00 I I ‘ I I ‘ I I 0 O ’ 0 I 6.00E+OO .O ‘ ‘ O O . I o to I I ‘ O 0 ‘000 I I 0 ‘ ‘ 0‘ O I ‘ 0— . 9 0- . ‘-0-‘ ”‘00000000000 -1.50E+OO . ‘ O ‘ ‘ ‘00000 . I 0 ‘ 0 O ‘ I I 0" ‘ O O I I ” 0 I I 0 0 I -9.00Ei-OO . ‘ ‘0 . I 0 ‘ 0 I I I I O I I“ O I -1.GSEO-01 . . . . . . DY '- 1.SOEO'00 UNITS PER LINE HIT (RETURN) TO (DN'TINUE... mm POST-PROCESSOR OPTION (G. T. P. E. X. (FULL)) :X PROCEED? (Y): mu READY TO LEAVE EUPORT6.0 ”‘ PROCEED? (Y): Euro a nice day. Appendix 2 Calling Tree for ENPOR16.0 This is s culling tree for ENPORTU.O us of August. 1984. It does not include IO utilities cannon to :11 routines. Procedures which sre shared with the linesr pucksge. ENPOR13.2. sre else not necessarily listed here. Those which hsve s note to see ENPORTB.2 sre those which are shsred wholly with the other pragrsn. MAIN EEDING GETERI INITLZ CL'S SETCOR INPROB UNSAVP INITLZ PROBDR GRAPDR ‘ GRPROC ‘ PARADR ‘ FUNCDR LSTFCN FCN'TR FNAXER FNSBl - FNSB23 SETFCN INIT DEFINE FNCDEF SOURCE ‘ SOLVDR DEOVL ICOND - SETIC PRINIC CENGIC A? TLIIG OUTVBL DISPDR LEAVDR FILEDR EELPDR NODSET Sec ENPOR15.2 ZOIC PRTOVL SETOVL SETALL NONLIN RKDRIV 44 (not currently implemented) GEIIN TPITXG TGET UGET SAVENL RUNGRT RESID UGET TPITXG TGET BPGET GETVAL Appendix 3 Subroutine Descriptions for ENPORTE.0 This is a list of short descriptions for the subroutines in ENPORTU.0 as of 8/09/84. This is not a production version. and is subject to change with virtually no notice. FUNCDR LSTFCN FCN'TR INIT DEFINE FCNDEF SOLVDR DEFOVL This is the driver for the description and definition of the constitutive and connective equations Prints the current definition of these functions. Prints a particular function definition. Controls the definition of a particular function. Initializes the nonlinear function description package. Describes the funcitonal dependancies in the system. Gets the user to describe the particular functions that he wants. This is the high level driver which controls selection of initial conditions. computational time limits. choice of display variables and problem integration. Selects some deault output variables by starting with the state variables. input variables and the derivatives of the state varaibles. It fills up the list with these and steps when full. 45 ICOND SETIC PRINIC CBNGIC ZOIC TLIR6 OUTVBL PRTVOVL SETOVL SETALL JSITX JSGET GETDN 46 Directs the setting. changing or displaying of the initial conditions. Allows the user to set all initial conditions. Prints the current initial conditions. Allows the user to set a particular initial condition. Sets all the initial conditions to zero. Gets the user to set the upper and lower computational time limits and the number of results saved or display. Directs the user to define the variables to be saved for later display. Displays the current varaible choices to be saved for later. Allows the user to choose a particular variable for later display. Gets the user to specify the whole list of things to save. Gives the user the choice of which integration method to use. The only active choice currently is a rather simple-minded Rnnge-Kutta method. Controls the computation of the time response using a simple third order Runge-Kutta technique. Sets up the static junction structure matrix. Eliminates the internal variables in the static junction structure matrix. Gets the variable names to be used to define the modulated functions. TGET UGET SAVENL RUNGKT RESID BPGET GETVAL 47 Sets up the modulated junction structure matrix. Eliminates the internal variables in the modulated junction structure matrix. Retrieves the current source field outputs. Saves the information needed for postprocessing at each time step. Uses a Runga-Kutta technique to find a state space trajectory. Describes the derivative of the state vector for the current time and state. Carves the junction structure matrix into the apporpriate subfields for later compuation. Returns a value for one of the nonlinear constituative “equations. Appendix 4 Source Listings for ENPORTE.0 C CSOLVDR ))))>>))))>)))>>>>>)) SOLVDR <<<(<<<(((<(<(<<<(<<<(<<<<<<(<((< C SUBROUTINE SOLVDR(PROCFG) C C-- SOLVDR CONTROLS THE SOLUTION MODULE FLO! C C-- PROGRAMMER: RCR'BERG SETTEIBER.1983 C iINsn'r wronncomNsLors momma SINsnr mmn'rmolmmmrs >INPI‘BK C CHARACTER‘70 STRING LOGICAL NEVLIN.ENDLIN.PROCFG.ISYES.READIT CHARACTER‘l ANS.DASH C DATA DASH /'-'l C c...smvnn..¢tOOOOOOOOO...O...O..00...OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO C C INITIALIZE XIN.TTN.TFIN...OUTPUT'VBL LIST. C CALL DEFINC(NX.XIN) CALL DEFTCP(MAXREX) CALL DEFOVL C IF (NX.E0.0) THEN WRITE(’.1000) 1000 FORMAT(/' "‘ NO STATE VARIABLES. NO SOLUTION.') RETURN ENDIF C C-- PRESENT THE SOLVER OPTIONS MENU lO CONTINUE CALL PAGE WRITE(‘. 1020) (DASH. I=1. 35). (DASH. 131.35) 1020 FORMAT(/' SOLVER OPTIONS'./1X. 35A1. 1 I' I: INITIAL CONDITIONS'. 2 /' T: TIME CONTROLS'. 3 /' S: SOLVE THE PROBLEM'. 4 I' F: SAVE THE RESULTS IN A FILE'. 5 /' H: HELP'. 6 /' X: EXIT nu: SOLVER ('DEFAULT) ' .IIX. 35AJ.) C 50 STRINGB' ENTER OPTION (X):' 48 49 CALL PROMPT(S'mING) NEILm-.TRUE. ANS-'X' CALL GET‘D(ANS.NEILIN.DIDLIN) C-- INTERPRET THE SELECTION AND ACT IT OUT IF (ANS.m.'I') THE! CALL ICOND(NX.XIN) ESEIF (ANS.m.'T') THE“ MAXREX - MAXRES DTSTR I MAXRES CALL TLIM6(TIN.TFIN.NSAV.D'TSTR) ESEIF (ANS.m.'S') mm CALL WTVBL CALL INTER ELSEIF (ANS.m.'F') THE! C CALL SAVRES ELSEIF (ANS.m.'H') mm IRITE(‘.'(/" FAT CHANCE FOR HH..P HERE...")') ELSEIF (ANS.m.'X') mm CALL PROCED(PROCFO) RETURN ESE 'RITE(‘.1090) 1090 FORMATU' INVALID OPTION SEECTD. TRY AGAIN.’ ) GOTO 50 EIDIF GOTO 10 C m C C CTLIM6 ))>)>>>)>>>>>)>>>>))>) TLIM6 <<<<<<<(<(<<(<<(<(((<(<(<(<(< C SUBWTINE TLIM6(TIN.TFIN.NSAV.DTSTR) C C-- EH6 SETS HE TIME coma. PARAMETERS FOR smUTION. C C-- VARIABLES: TIN INITIAL TIME C TFIN FINAL TIME C NSAV TOTAL NUBER STAGES TO STORE (INCL AT TIN) C DTSTR TIME INTERVAL BETWEEN! STORES C C WICAL NE'LIN.ENDLIN CHARACTER‘7O STRING C COOOTLIIGOOOOOOOOOO0.000000000000000000...0“...OOOOOOOOOOOOOOOOOOOOOOO C WRITE(‘.1000) 1000 FORMATU ' TIME CONTROL PARAMETERS FOR SOLUTION: ') C VRITE(‘.'(1X)') 50 'RITE(STRING.1010) TIN 1010 FORMAT(' mm INITIAL TIME ('.1PEIZ.4.') :') CALL mourmsnmo) NEILIN=.TRUE. TLOIO. mI=1.E10 CALL GEIRL(TIN.TLO.'mI.NE'LIN. mum) 'RITE(S'TRING.1020) TFIN 1020 FORMAT(' mm FINAL TIME ('.1PE12.4.') :') CALL mournsnmo) NE'Lm'.'mUE. CALL GEIRL(TFIN.TIN.THI.NEILIN.E(DLIN) IRITE(STIING.1030) NSAV 1030 FORMAT(' EITER NUIBFR OF POINTS TO STORE ('.I3.') :') CALL PROMPT(STRING) NEVLIN-JRUE. ILO=2 IHIBDTS'm CALL GETIN(NSAV. ILO. IHI.NEILIN. ENDLIN) C DTSTI'l (TFIN-TIN) / (NSAV-1) C C WRITE(‘.1040) DTSTR C1040 FORMATU ' THE STORAGE INTERVAL IS: ' .1PE12.4) C RETURN ETD C C CINTEG >>>>>>>>>>>>>>>>>>>>>>>> INTEG <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< c sonaooTINB mm 2— mm: comms mnuoN nou'rINE SELECTION glem'r wroumouoNnLoxsmnocsx SINSERT mmar>coamNsLoxs>smmmc C LOGICAL NEILIN.ENDLIN CHARACTER‘70 STRING CHARACTER‘l ANS.DASH C DATA DASH/'-'/ C ceeeINTEBeeeeeeeeeeeeeeeeseeeeeeeeeeeeeeeeeeeeeeeseeeeeeeeeeeeeeeeeeeee C C-- PRESENT THE.INTEGRATOR MENU 10‘ CONTINUE lRITE(‘.1020) (DASH. 131,45) . (DASH. I81 .45) 1020 FORMATU' INTEGRATOR OPTIONS' ./1X.45A1. 1 / ' A: ADAMS-BASHFOURTH PREDICTOR-CORRECTOR INTERATION ' . 1090 C C MM§N 51 : BACI'ARDS-DIFFEREV'TIATION METHOD' . : OTHER EXOTIC METHODS’ s RUME-KUT'TA (THIRD-ORDER) ' . : RETURN (IDEFAULT)'./1X.45A1) /' ’0 ll ’0 DEOU SHINGI' m OPTION (Q):' CALL PROMPT(STRING) NEILINIIJTRUE. ANS-'Q' CALL GEI'D(ANS.NE'LIN.EVDLIN) IN'IBBPBB'I nus (Dom: IP (ANs.Bo.'A'> THEN CALL NmLIN('A') mm as ' PRINT -. ' 'mIs man 15 NM (mm: AVAILABLE. ' PRINT an ' ELSEIF (ANs.Bo.'B') 1mm CALL NCNLINUB') PRINT tn ' PRINT t. ' 1313 won 13 nor cum: AVAILABLE.’ PRINT on ' n.3BIP(ANs.Bn.'B') THEN CALL RKDRIV ELSEIF (ANS.EO.'O') “mm PRINT on ' * mmvuu" ARE you KIDDING?")') PRINT tn ' ESEIF (AN8.Bo.'o') 1mm 2mm use WRITE(*.1090) PoBNA'I'u' INVALID OPTION. PLEASE n! AGAIN.') mDIP STRING!” mm OPTION OR HIT (BET): ' CALL PROMPT(STRING) NE'LIN-.'TRUE. ANS-0 I CALL GET'D(ANS.NE'LIN.E{DLIN) IF (ANS.EO.' ') THE! 6011) 10 ESE GO'IO 50 DIDIF HUD CENDINTEG C CRKDRIV >)>>>))))>>>>>>>>>>>>>>>> RKDRIV <<<<<<<<<<<(<(<<((<(((<<<<(<<( C SUBRWTINE RKDRIV 52 DOUBLE PRECISION Y(10) .DY(10) .SLOPE(10) .DPTIM WEB LBS(12) COMI'JN llBS/IBS . LINDIS. LINSTO. LINJ'S LCXiICAL LINDIS. LINSTO. LINJ'S. NEWLIN. ETDLIN. SUCESS. TRACE. NE'FCN CHARACTER‘72 STRING SINSERT FNPOR'DCONJDNBLOKS >UTILBK SINSBB'I I-NPORT)COIODNBLOKS>SYBGBK iINSBB'r DIPORT>COMDKJNBLOKS >CLASBK SINSERT I-NPORT)COIOI)NBLOKS >BDUCBx $INSBB1' INFORT>COMLDNBLOKS>SOLNBL6 SINSBB'I mmRT>COMIDNBLOKS >COIIAREA C 0000 OOGG 000 000000 NE'LIN II .‘mUE. T'RACE - .FALSE. Set up LBS array of pointers in the S array... CALL J’SM'TX PRINT ‘.' IERRF 8'.IERRF IF(IERRF .NE. O)RETURN CALL JSGET(NBEX.WIN..FALSE.) PRINT ‘.' IERRFs'JERRF IF(IERRF .NE. O)RETURN CALL PAGE CALL PROIIP'l‘(' Want a linear dissapation field? (YES)') LINDIS . .FALSE. CALL YORN(LINDIS) CALL PROMPT(' Want a linear storage field? (YES)') LINSTO I .FALSE. CALL YORN(LINSTO) CALL PROMPT(' I suppose you"re happy with linear two-ports?') call prompt(' (YES)') LINJ'S I .FALSE. CALL YORN(LINJS) Here we define the input to the funciton which determines the modulus of the modulated transformers and gyrators. IF(.NOT. LINJ'S)THm DO 9 I=l.NFT PRINT ‘.' For two-port number'.I.'. which state variable' PRINT ‘.' do you want to use to determine the modulus? (1)' MDINU) t 1 CALL GETIN(TIDDIN(I).1.NFT.NEWLIN.MLIN) 00000000060 53 TMTYPE(I) I 'X' CONTINUE ENDIF CALL TPMTX6(SUCESS.LINJS) CALL TGET(NBEX.NFT.TRACE) Now the dissapation field.. CALL FLMTX And the storage field.. CALL FSMTX(IFSZ.IFS3.IFS4) lRITE(STRING.'(" Ihat time increment do you want to use". 1 ". for computation? (".1PE9.2.")")')DTSTR CALL PROMPT(STRING) DELTA I DTSTR CALL GETRL(DELTA.1E-15.DTSTR.NEWLIN.ENDLIN) DT I DELTA TIME I TIN DPTIME I TIN STOTIM I TIN + DTSTR CALL UGETCTIME) DO 10 I I 1.NFI Y(I) I XIN(I) CONTINUE I I 1 CALL RESID(NFI.DPTIM.Y.SLOPE.DY.IERR) CALL SAVENL(I.Y.DY) CALL CLOCX(START) CALL BDNCB'I(NPI.TIIB.Y.DY.DI) PRINT a! rue-um IP(TINB.Bn.smm)TBm I - I + 1 3mm - smm + DTSTR DP'I‘IM a mu: CALL BPBIMNPI, DP'I'IN.Y. SLOPE.DY. IBBB) CALL SAVM(I.Y.DY) DT . DELTA mm TIME I TIME + DELTA C 54 IF(TIME.GT.STOTIM)THEN DT I TIME - STOTIM TIME I STOTIM ENDIF IF(I.LT.NSAV)GOTO 20 CALL CLOCX(FINISH) PRINT ‘.FINISH-START. 'CPU seconds used.’ RETURN END CENDRRDRIV )>>>)))>)>)))))))))))>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C CRUNGET )))>>>))>))>)>>))>>>>))> RUNGKT (((<<<<(<(<<<<<(((<(<<(((<<<<< C C 50 SUBROUTINE RUNGRT(N.TIME.Y.DY.DELTA) REAL TIME.DELTA.XO(10).K1(10).K2(10) DOUBLE PRECISION T.Y(‘).DY(’).SLOPE(10).R(10) DATA SLOPE /10'0.0D0/ T I TIME CALL mm>>)>>))>>>>>>>>>>>>>>>>)>>)))>>>)>))))>>>>>>)>>>)>>)>)>>>> C CSAVM >>>>>>>>>>>>>>>>))>>>>>> SAVEU. ((<<<(<(<<<(((<<<<((<<<<<<<<<<< C SUBRWTINE SAVM(NS.Y.DY) c C-- SAVPNL: SAVEO x.U.Dx AT STAGE Ns IN RES(NS.NV)FOR NONLmEAB VERSION c C M:- STAGE FOR smBINO RESULTS c Nv- NUIBER OF VBL m OUTPUT LIST (IN ORDER x.U.DX) C DOUBLE PRECISION Y(‘).D’f(‘) 31mm mmaT>COENONBLors>smNBx.6 c COCOSAVN‘COO0......OOOOOO0.0QOOOOOOOOOOOOOOOOOOOO0.0000COOOOOOOOOOCOOO C Nv-I C no 20 I=I,Nx IF (XSAVL(I)) TEEN RES(NS.NV)-Y(I) Nv=Nv+1 1P (NV.OT.NBEO) BETUBN ENDIF 20 CONTINUE C D0 40 I-I,NU IF (USAVL(I)) mm mmsmw-Um Nv-Nv+1 , IF (NV.OT.NBEO) BETUBN DIDIF 40 CONTINUE C no so II1.N! IF (DSAVL(I)) Tam RES(NS.NV)IDY(I) NV-NV+1 IF (Nv.cT.NBEO) RETURN ENDIP . so CONTINUE C BETUBN END C C SUBROUTINE RESID(N. TIME.Y. SLOPE. R. IRES) 56 DWBLE PRECISION R(N).EOPE(N).Y(N).TIME REAL ZI(10).DI(10).DO(10) INTEGER IPZI(10).IPDI(10).UBS(12).IPDO(10).LENDI.DPOINT COMIDN ILBS/IBS .LINDIS.LINSTO.LINJ’S LOGICAL LINDIS.LINSTO.LJNJS.TRACE.SUCCES SINBEBT ENPOBT>CONIONBLoxs >SYBGBK iINSEBT ENPOBT>COEEONBLOKS >CLASBK SINSEBT ENPOBT>CONnnNBLors>aDUCBx SINSEET mPOBT>coumNBLoxs>smmr SINSERT mPOBT>CONNONBLoxs>CONABEA C OOOOOOOOOOOOOOOOOOOOOOOGOOO INTEGER IPU(MAXU).IPX(MAXX).LENX DATA TRACE /.FALSE./ Here we define the implicit differention equations that make this whole program worthwhile. The implicit formulation which must equal zero to solve the problem is set up and the 'residule' which is left is given to integration package to chew on. In general we will be doing the following steps: 21 I FSll‘XI + FSTZ‘XD DI I 831 ‘21 + 834 ‘U (explicit case) DO I FL'DI (linear case) or DO I PEIL(DI) (nonlinear case) R I Sll‘ZI + 813‘DO + Sl4‘U - AIS Ihere R is a residual that LSODI forces to zero. and SLOPE is the current approximation of the derivatives of the state variables (both independent and dependant). In the case of the dependant state variables there are also algabraic constraints with must be tacked on to the bottom of the vector R. First we have to define the non-linear two port elements (MTF and MGY) 06060 and reduce them into the junctin structure matrix (8). Go get the inputs... SHORTITIME CALL UGET(SHORT) DO 3 II1.NU IPU(I)II CONTINUE LENUINU 57 Stuff current state into XI and XD... GO LEIX I NFI DO 5 II1.NFI X(I)IY(I) IPX(I)II 5 WE IF(.NOT. LINJS)THEN CALL T'Pfl'X6(SUCCES.LINJS) IF(.NOT. SUCCES) THE] PRINT ‘.' TPM'TX6 NOT SUCCESSFULL. CALLED FROM RESID... DIDIF ENDIF CALL TOET(1‘8EX.NFT.TRACE) CALL BPGET LBS(1)I1 DO 10 MI1.11 0 LBS(M+1)IIBS(M)+LS(M) First. F811‘XI OOOH IF(LINSTOYTHDI 0 CALL CSMPY (NFI.NFI.LE{FS. IPFS. FS.1 .LEIX. IPX.X) CALL GEIVX(ZI.LETZI. IPZI) on ESE 21 I PHIS(X) OOOO DPOINTIT LMIINFI D0 72 1.11m!) IF(ABS(IBT(I)).m.1 .OR. ABS(IBT(I)).EQ.2)THW D0 330 .7131me IF(IEQZS(JJ.1) .m. I) Tum NOMIJJ 60m 633 MIR 330 CONTINUE 633 1121 (DPOINT)IDPOINT CALL GETVAL(X(IPX(DPOINT) ) .ZI (DPOINT) .NE. 1 IFCNTP(NOEQ)) DPOINTIDPOINT-I-l ENDIF 72 CONTINUE ENDIF C C Now SSl‘ZI 000 00606 000 GOG COO 0006 OOOGOO 58 CALL CSMPY(NFL.NFI.LS(9) .IPS(IBS(9) ) . S(LBS(9) ) .1.LMI. IPZI.ZI) CALL GETIX(DI.LENDI.IPDI) $34 ‘ U... CALL CSMPY(NFL.NFS.LS(12).IPS(LBS(12)).S(UBS(12)).1.LENU.IPU.U) And now add them into DI.. put S34‘U into D0 (temporarily)... CALL GEI‘IK(DO.LE{DO. IPDO) And now add them.. CALL CSADD(LENDI.IPDI.DI.LENDO.IPDO.DO) Put the result in DI.. CALL GET'X(DI.LENDI.IPDI) Here we calculate the DO vector as a function the DI vector for either the linear or nonlinear explicit case. IF (LINDIS mum DO I FL'DO CALL CSMPY(NFL.NEL.LENFL.IPFL.FL.1.LENDI.IPDI.DI) CALL GETIK(DO.LENDO.IPDO) ELSE DO I PEIL(DI) THE FOLLOWING IS A KLUGE TO GET PAST'THE INABILITY TO REFERENCE A.ZERO ELEMENT OF A SPARCE MATRIX... DO 40 KI1.NFL DO(X) I O. CONTINUE DO 50 XI1.LENDI DO(IPDI(X)) I DI(X) CONTINUE DO 660 XI1.NEL DI(X) I DO(X) IPDI(X) I K 660 CONTINUE DPOINTIl LENDOINFL 30 60 GOO GOO GOO 000 59 DO 7 TILED IF(ABS(IBTTI)).EO.3)THEN DO 30 JJITJBD IF(IEQZS(JJ'.1) .m. I) THE! NOMIJJ' GOTO 60 MDIF CONTINUE IPDO (DPOINT)IDPOINT CALL GETVAL(DI (IPDI (DPOINT) ) .DO(DPOINT) .NOm. 1 IFCNTP(NOEQ) ) DPOINT‘IDPOINT+1 MIF CONTINUE NIP 21 I Sll‘ZI (temporary storage in 21) CALL CSMPY(NFI.NFI.LS(1) .IPS(LBS(1) ) .S(TBS(1) ) .1.LDIZI. IPZI.ZI) CALL GET'X(ZI.LE(ZI.IPZI) DO I Sl3‘D0 (temporary storage in D0) CALL CSMPY(NFI.NFL.LS(3) .IPS(LBS(3) ) .S(LBS(3) ) .1.LmDO. IPDO.DO) CALL GET‘X(DO. LEWDO. IPDO) DI I Sl4‘U (temporary storage in DI) CALL CSMPY(NFI.NFS.LS(4) .IPS(TBS(4) ) .S(I.BS(4) ) .1.LmU. IPU.U) CALL GET'R(DI.LEWDI. IPDI) Here we add the various components of XDOT together... CALL CSADD(LE(ZI. IPZI. ZI. WDO. IPDO. DO) CALL GET'K(ZI.LMI. IPZI) CALL CSADD(LMI. IPZI. ZI. LEUDI. IPDI.DI) DO 11 IIl .N R(I)I-EOPE(I) CONTINUE DO 20 II1.N B(M(I))-u(1)+B(wa(I) ) CONTINUE RETURN ETD SUBRWTINE 1 AC RETURN m 6O (5 SUBMTINE GE'NULOCLJEOEJPLOCL) This subrououtine takes the contents of the work array VI and puts it in the array LOCI. along with the pointer array IPLOCL and the length count LNLOCL. INSERT EWPORT>COMUNBLOXS )RDUCBX . fiflifififififl REAL LOCL(‘) INTEGER IPLOCL(‘) DO 10 II1.N IPLOCL(I)II'X(I) LOCL(I)I'X(I) 10 CONTINUE LEOE I mex RETURN BID CDCDCDCS CFUNCDR >>>>>>>>>>>>>>>>>>>>>>>> FUNCDR <<<<<<<<<<<<<<<<<<<<<<<<<<< c SUBROUTINE FUNCDR(PROCFG) c c— PBOOBANIUIB: c. ALLEN JUNE 1984 C C-- PURPOSE: FUNCDR' OONTBOLS TEE FUNCTION SPECIFICATION OPTIONS C sINSEBT mPORT>COMlnNBLOKS>SYBGBK SINSERT ENPOBT>coumNBLoxS )UTILBK SINSEBT ENPORT>COMIDNBLOKS >EDUCBx SINSERT mPOBT>CONmNBLOKS>CONABEA C GARACTER‘7O STRING LOGICAL NE'LIN.E(DLIN.PROCFG.FULL CHARACTE‘l ANS.DASH CHARACTER‘SZ FNAME C DATA DASH/ '-'/ C COOOFUNCDROOOOOOOOOOOOOOO...00......QC.........OOOOOOOOCOOOOO0.00000000.0 C IF (NE.m.O) THEN! IRITE(‘.1000) 1000 FORMATU' SYSTEM GRAPH EMPTY. NO ACTION. ') PROCFGI.FALSE. RETURN ENDIF FULLI.TRUE. c... 10 1020 mu‘uNH 61 PRESM IRE FUNCTION OPTIONS mu CWTINUE WRITE(’.'(1X/1X)') IF (FULL) THEN 'RITE(‘.1020) (DASR.I-1. 41). (DASH. I-1.41) FORlAT(/' FUNCTION OPTIONS'./1X. 41A1. /' L: LIST CURRENT FUNCTIONS'. /' C: CHANGE SEECTW FUNCTIONS'. /' S: SET ALL FUNCTIONS'. /' U: UNIFORI FUNCTIONS (ZERO OR (NE) 'o /' R: RELP'. I ' X: EXIT FUNCTIONS SECTION ('DEFAULT) ' ./1X41A1) STRING" m OPTION (X): ' ELSE STRING" m FUNCTION ORION (Lflclslu'n010m>):' MIF READ 13E SELECTION CALL PRONPNS‘RING) NE'LIN-JRUE. wk! 0 CALL GET'DMNS. NE'LIN.mDLlN) INTERPRET THE SELECTION AND ACT IF (ANS.m.' ') THE] IF (FULL) mm ANS-'X' ELSE FULL-(RUE. GOTO 10 DIDIF WDIF IF (ANS.m.'L') THEN CALL LSTFCN CALL GOCN ESEIF (ANS.m.'C') THEN 'RITE(‘.'(/." FEAT FUNCTION DO YOU WANT TO CRANGE7")') READ(‘.‘)NOF CALL SE'I'FCN(NOF) ELSEIF (ANS.m.'S') THEN CALL INIT IPNTU) ' 1 CALL DEFINE OPTION '3 'U' CALL FCNDEF PRINT ‘,' ' CALL GOG! ELSEIF (ANS.m. 'U') mm PRINT ’.' NUT CURRENTLY AVAILABLE. SORRY.’ PRINT ‘, ' ' ESEIF (mam. 'R') THEN FNANE- ' PUBLIC >85 .2 )REP. FUNCTIONS ' 62 CALL RRFILE(FNANE) ESEIF (ANS. m.'X') TEE] CALL PROCED(PROCFG) 11mm ELSE 'RITE(‘.1090) 1090 FORIATU' INVALID OPTION SEECTED. IR! AGAIN.') MIF FULL-.FALSE. GOTO 10 C m C SUBRGTTINE LS'TFCN C Cwm B! : C. m JULY 1984 C C 11113 ROUTINE LIST ALL OF TEE cum: DEFINED FUNCTINS C . SINSEET ENFOET>CONIONELOTS >EOUCUK SINSEET EUPORT>CORDNBLOKS>COIIAREA c c CALL PAGE C IF(NES.m.O)TRm PRINT ‘. ' NO NNSUTUTNE mUATIONS CUM! DEFINED. ' PRINT " RETURN DIDIF D0 10 K '3 1.Nms IF(W'HNK) .LE. 6) mm IF(EVS(K) .BQ.0)'mm 'EITE(‘;9020)K. m10(x.1).NFCN(K).FQNSE’1‘(IFCN’I'P(KH. ELSE IF(RV8(K).m.1)THEN IRITE(‘.9020)K. GIO(K.2).NFCN(K).FCNSET(IFW(K))o + GIO(K.1).NFCN(K) m1]? IE HAVE A NUL'TIPORT. .. GOO ESE IF(RVS(X) .m.0) TEN 'RITE(‘.8820)X. CUIO(X.1) . IE028(K.1) .FCNSET(IFCNTP(X) ) . + GIO(X.2).IEQZS(R,2). + CRIO(K.3).IEQZS(X.2).FCNSET(IFCNTP(X)). + (BIO(K.4).IEQZS(X.1) ESE 'RITE(‘.8820)X. CRIO(K.2) . IWS(K.2) .FCNSETUFCNTPK) ) . + GIO(R.1).IEQZS(K.1). 63 + CRIO(X.4).IEOZS(K.1).FCNSET(IFCNTP(X)). + GIO(X.3).IEOZS(X.2) C MIF C MIF 10 CONTINUE C C—--FORNAT-----—- C 9020 FORMAT(5X,I3.': '.A1.IZ,' - '.A6.' ('.A1,12.')') 8820 FORNA'T(5X,I3,': ',A1,IZ,' ' '.A6.' ('.Al.12.')'. + I.SX.A1.IZ.' - '.A6.' ('.A1.12.')') C RETURN HUD C C C C SETFCN )>>)>)>)>>)>))>>)))))))>>)>>))))>>)>))))>>)))>))>>>>>>>)>>)>>)> C SUBRWTINE SETFCN(NOF) C C IBIS RWTINE ALLOIS THE USER 1U CHOOSE 'IRE FUNCTIONAL DEPENDWCY C FOR THE CONSTITUTIVE MPS'TRON AS DETERNINED BY NOF. C C mm B! G. ALLEN JUL! 1984 C LmICAL WE. REP.NE'LIN.E{DLIN INTWER POINT. CIPS'H. C CHARACTER ‘72 STRING CRARACTER‘M ANS C SINSERT INFOET>CONEONELOTS >RDUCBK SINSERT INPOET>CONNONELOIS>CONAEEA C WE - .FALSE. NE'LIN '- .IRUE. C 'RITE(‘,'(/." FUNCTION",IZ," NW DEFIN- AS ")') NOF CALL FCN'n(NOF) C mm(.'l(ll 00)!) 10 WRITE IF(ANS .m. FCNSET(IFCNTP(NOF)))TREN RETURN WDIF 64 C POINT ' CNPS'IR(ANS.FCNSET.IAXFUN) C IF(POINT .m. 0)'mm 'RITE(‘.'(//.A." IS NOT A FUNCTION TYPE EAT IS (EMILY ". 1 "DEFINED. ".//)')ANS 'RITE(S'mING. '(" WELD YOU LIKE TO SEE A LIST? (NO)")') CALL PROIPT (STRING) CALL YORN(RELP) C IF(REPYTREQ PRINT ‘ PRINT ‘. ' TRIS IS A LIST OF FUNCTIONS. . . ' ENDIF ESE C C THE USER GAVE US AN ACCEPTABLE FUNCTION TYPE.. . C CALL FNAKER(NOF.POINT) CALL WNNOF) C DONE - .TRUE. C MIF C IF (.NOT. WE) GOTO 10 C RETURN BID C msmcrOOOOOOOOOOOOOOOOOCOO...0.000000000000000.00000000000000000000C C CCNPS'R)>))>>))>)>>>>>>>>>>>)>>>>)>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>)C TRIS FUNCTION RETURN A POINTER TO THE CHARACTER ARRAY ’LIST' WHICH INDICATES TUE EENENT INEE ARRAY HAT NATCRES THE (BARACTER VARIABLE 'EXANPL'. (NLY THE FIRST FOUR (IARACTERS ARE CONPARE, AND A ZERO IS RETURN. IF NO HATCH IS FOUND. mum BY G. W JULY 1984 (5000005000 INTEGER FUNCTION CIPS'TR(EXANPL. LIST. LEON) fl CRARACTER‘(‘) 'EXANPL.LIST(‘) CIPSTR - 0 DO 10 1'1..me IF(EXANPL .m. LIST(I)(1:4)) TERI CMPSTR - I RETURN ENDIF 65 CENDCIPSTR)>>)>>)>>)>>>>>>>>>>>>)>>>>)>>>>>>>>>>>>>>>>>>>>>>>)))>>>)>>>C 10 CONTINUE RETURN END C C C--FIAXER C C C---'RITTEN BY : C SUBROUTINE.FIAXER(NOF.FCNTYP) G.ALLEN JULY 1982 SINSER‘I‘ mmnnoouomwxmnwmx Smsm'r mmumomunwxs mom C OOOOOOOGOOGOOOOOOOOOOOOOOOOOOOOOOOOOn INTEGER FCNTYP VARIABLE : IFCNTP -- STORE FUNCTION TYPE BELOI. ‘9”. ALL FUNCTION INDICES USED IN OEER SUBROUTINES ARE ”“9 ACCORDING TO THE FOLLOWING TYPE NUDE. TYPE DEFINITION 01 CONST Y-CO 02 LINEAR Y-CO+C1‘X 03 SIN Y-Cl‘SIN(C2‘X+C3) + C4 04 COS Y'Cl‘COS(C2‘X+C3) + C4 05 ASIN Y=C1‘ASIN(C2‘X+C3) + C4 06 ACOS Y'Cl‘ACOS(C2‘X+C3) + C4 07 ABSZ Y-CI‘X‘ABS(X) 08 ABSQRT Y-Cl‘SIGN(X)‘SQBTKABS(X)) 09 EXP Y'CI‘EXP(C2‘X+C3) + CH 10 ALOG YtCl‘ALOG(C2‘X+C3) + C4 11 DSIN Y-Cl/SIN(C2‘X&C3) + C4 12 DCOS Y=C1/COS(C2*X+C3) + C4 13 DASIN Y-C1/ASIN(C2’X+C3) + C4 14 DACOS Y-CI/ACOS(C2‘X+C3) + C4 15 POLY Y-CO+C1‘X+C2‘X“2+C3‘X“3+C4‘X“4 16 COLUIB Y=C1‘ SIGN(X) 17 BRKP'T Y-Y1+C1‘ (X-Xl) . X (.Xl Y-Y1+C2‘(X-X1) . X1< X 18 SATUR Y=Y1 . X <‘ X1 Y=Y1+((Y2-Y1)/(X2-X1))‘(X-X1) . X1