THESIS

 

This is to certify that the

thesis entitled

A PARALLEL COMPUTATION METHOD
FOR DYNAMIC SYSTEMS
WITH COUPLED NONLINEAR DISSIPATION

presented by

Tong Zhou

has been accepted towards fulfillment
of the requirements for

Masters (hxgeein Science

 

 

,1 ”W. C
- " k,'
(K4 '

Major professor

 

Date fl/Laéfl WW?

0-7639 MSU is an Affirmative Action/Equal Opportunity Institution

 

 

 

MSU

LIBRARIES
——.

‘r

 

 

RETURNING MATERIALS:
Place in book drop to
remove this checkout from
your record. FINES will
be charged if book is
returned after the date
stamped below.

 

 

 

 

 

A PARALLEL COMPUTATION IETEOD
FOR DYNAMIC SISTERS
WITH COUPLED NONLINEAR DISSIPATION

By

Tong Zhou

A THESIS

Submitted to
Iiehigan State University
in partial fulfillment of the requirements
for the degree of

IASTER OF SCIENCE

Department of Mechanical Engineering

1984

ABSTRACT

A PARALLEL COMPUTATION IETEOD
FOR DYNAMIC SISTERS

'ITE COUPLED NONLINEAR DISSIPATION

Tbng Zhou

Nonlinear algebraic loops in system equations may prevent the
subsequent reduction of the equations to an explicit state-space form.
It may make system simulation very difficult to accomplish. The inci-
dence of nonlinear algebraic loops in mathematical system equations
often is associated with the existence of nonlinear dissipative

effects in physical systems.

A computation method for nonlinear algebraic loops stressing par-
allelism has been developed. The bond graph augmentation method
converts an algebraic loop field into a dynamic subsystem that exhi-
bits the proper static characteristics at steady state and employs a
two-time-scale integration technique. In seeking computational effi-
ciency. minimizing the augmentation order and Optimizing the parameter
selection play key roles. An augmentation sequence and a general rule
for parameter selection for arbitrary n-th-order subsystems have been

suggested and numerically tested in several cases.

To
my motherland China.
my parents I. J. Zhou . B. X. Lee.
and

my wife Shuling

ACKNOWLEDGEMENTS

The author vishes to express sincere thanks to Dr. Ronald C.
Rosenberg for his guidence and encouragement over the past years, and

for serving as major professor.

Appreciation is also expressed to the Albert R. Case center for

providing computer facilities for this research.

CONTENTS

LIST OF TABLES
LIST OF FIGURES
1 INTRODUCTION
1.1 The Problem of Algebraic LOOps in System Equations ..........1
1.2 Previous Research 'ork ......................................2
1.3 Bond Graphs and Algebraic LOOps .............................6
1.4 An Approach Stressing Parallelism ..........................13
2 BOND GRAPH DYNAMIC ADGNENTATTON NETEOD
2.1 Identification and Partitioning of R-fields ................15
2.2 Dynamic Augmentation of Subfields ..........................18
2.3 Tho-time-scale Integration .................................34
3 SOME COMPUTATIONAL CONSIDERATIONS
3.1 Effect of Eigenvalues on Computational Efficiency ..........37
3.2 Parameter Selection in Linear erields .....................41
3.3 Parameter Selection in Nonlinear R-fields ..................62
4 NUMERICAL EXAMPLES
4.1 Linear Case ................ ..... ...........................70
4.2 Nonlinear Case .............................................83
5 CONCLUSION
5.1 Summary ....................................................90
5.2 Future Development .........................................92
REFERENCES ......................................................94
APPENDIX A: Program Calling Tree ................................96

APPENDIX B: Program Listings ....................................97

LI ST OF TABLES

TABLE 3-1 EJGENVALUES
TABLE.3-2 EJGENVALUES
TABLE 3-3 EIGENVALUES

TABLE 3-4 EIGENVALUES

LIST OF TABLES

OF A SECOND-ORDER AUGNENTATION...............48
OF A THIRD-ORDER AUGNENTATION................56
OFMMIMY SUBSYSMQQO00000000000000.0058

OF A SUBSYSTEM WITH CHAIN STRUCTURE .......61

LI ST OF FIGURES

Figure
Figure
Figure
Figure
Figure
Figure
Figure
Figure
Figure
Figure
Figure
Figure
Figure
Figure
Figure
Figure
Figure
Figuer
Figure
Figure
Figure
Figure

Figure

LIST OF FIGURES

Bond Graph Structure and Key Vectors ...................8
Refield in a Physical System ..........................10
Physycal Systems with a Coupled R-field ...............16
The Partitioned Subsystems ............................17
The Interactions Among the Partitioned Subsystems .....19
A Static Field and its Augmentation ...................21
Schematic Diagrams ....................................23
Partitioning ..........................................26
Causality Orientations ................................29
An Example ............................................32
Schematic Diagram of Two-time-scale Inteegration ......36
Diffusion Convergence ...............................39
First-order Augmentation ..............................42
Second-order Augmentation .............................42
The Eigenvalue DIstribution ...........................50
Chain Junction Structure and its Augmentation .........53
An R-field : Example 1 ................................55
An R-field : Example 2 ................................57
A Third-order Augmentation ............................60
A Nonlinear R-field ...................................64
The Global Bond Graph Nodal for the Numerical Example .71
The Augmentation Result for the Example ...............72

The sub‘y‘ton' OO...O0.0.0.0.0...00.0.0000000000000000073

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

4-8

4-9

4-10

4-11

4-12

4-13

Local-time Integration over the First Global-time step 77
Local-time Integration over the Second Global-time

Step ..................................................78
Numerical Result of Linear Case Obtained by the
Two-time-scale Integration ............................79
Dynamic Response Obtained by the Two—time-scale
Integration ...........................................80
Numerical Result Obtained by the ENPORT-S .............81
Dynamic Response Obtained by the ENFORT-S .............82
Local Integration with Arbitrary Parameters ...........86
Local Integration with Optimal Parameters .............87
Numerical Result Obtained by the Two-time-scale
Integration ...........................................88

Dynamic Response of the Nonlinear Case ................89

1 INTRODUCTION
1.1 The Problem of Algebraic Loops in System Equations

A dynamic system consisting of a number of lumped-parameter ele-
ments may be described by ordinary differential equations in which
time is the independent variable. By using vector notation. an n-th
order differential system mmy be expressed as a set of first-order
vector differential equations. This is the state-space representa-

tion.

Because the modern trend in engineering systems is toward greater
complexity due mainly to the requirements of performing complex tasks
with good accuracy. the state-space representation has come to play a
very important role in modern system control and system dynamics. The

general form of explicit state-space equations is

in) = GF( xmmm ) (1.1)

where the X is an n-dimensional vector, U is an m-dimensional vector,
and G is a vector function. The superdot denotes time differentia-

tion.

If the system is linear. Eq.(1.1) takes on a simpler form, as

shown below

i = [AM + [810 (1.2)

where x and U are defined as before and [A] and [B] are matrices of

appropriate dimensions.

In formulating equations for a dynamic system the explicit
state-space equations often are preferred. If there are intermediate
variables which are introduced in the initial problem formulation,
there is the possibility that algebraic lOOps may exist in the system
of equations. This situation may make it impossible to find an expli-
cit form. More generally. the system equations can be written in the

form

in) = apt x. u. n ) (1.3.)

R(t) - 61} x. u. a ) (1.3s)

where B is a vector of dimension r.

For a linear constant-coefficient system. an explicit state-space
representation usually can be obtained by eliminating the intermediate
variables from the equations. However. if even one of the G functions
is a nonlinear function. it might make the elimination of B difficult.
or even impossible. The existence of nonlinear algebraic loops in
system equations may prevent subsequent reduction of the system equa-
tion set to an explicit state-space form. That leads us to search for

efficient computational techniques to treat such problems.

1-2 Previous Research 'ork

The problem of algebraic leaps has been studied by many research-
ers and many simulation programs are available that will diagnose the
existence of algebraic loops in the equation set. Among those are
CSNP III. CSSL-IV. DARE. and SCEPTRE [1.2.3.4]. Operationally. a loop
diagnostic occurs following the equation sorting process. In this
process.as the system equation set is manipulated . mutual algebraic
loops are identified. Typically. execution of the program is termi-

nated and appropriate modifications must be performed.

lost methods of numerical solution for nonlinear algebraic
equations are extended from those used for solving linear systems.
The Jacobi method. the Gauss-Seidel method . the successive overrelax-
ation method and the symmetric successive overrelaxation method are

exampleslSJ.

Among the methods for linear systems one which may be applicable
to the nonlinear case is the use of an iterative procedure. each step
of which involves the solution of linear algebraic equations. For

example. consider the system

4n, - n, + (1 I10):“1 = 1 (1-4')

'u, + 4u3 + (1/8lu: = O (1.4b)

An iterative method can be designed such that

)
4uin+1) ' u§“+" ‘|’(l/10)e“(n B 1 (1.5a)

-u§n+" - 4u§n*" +(1/s)<u‘n’)’ = o (1.5s)

T0 determine “:n+1) and n§‘*1’ from previous uin) and ugn)

involves solving a system of two linear equations with two unknowns.

Another method is the Newton method for solving a nonlinear sys-

tem. Consider the case

f;(u; .u, ) 8 0 (1.6a)

f.<n. .u. > . 0 (1.6b)

which we may write as Fu = 0.

The Newton method is defined by

“(n+1) = n‘h) _ (F'u(n))Fu(n) (1.7)

where F'u‘n) is the Jacobian matrix.

The Jacobian method is one of the basic methods for conducting

the iteration. For a nonlinear system of equations

f1(n1 .113 an, ) g 0 (1e8.)
f.(n. .u. .u. ) = 0 (1.8b)
fl(u1 an: an! ) g 0 (1.80)

.n:n+1) 'n:n+1) by

We find u£n+‘)

f1(u£n+1) .ugn) .ugn) ) = 0 (1.9a)
£,(n§n’ .u§n*" .ufin’ ) = o (1.9s)

f,(u:n) .ugn) .u£n+:) ) = 0 (1.9c)

At each time step . one has to solve a single nonlinear equation

for one unknown.

The Gauss-Seidel method is the same as the Jacobi method except
that at each time stage one uses the latest available values. Thus

for the above case

f.(u§“*" .ufin’ .nfin’ ) - o (1.10.)
f.<n§“*" .n§“+" .uin’ ) - o (1.10s)
f.(n§“+’) .n§‘*" .nS“*" ) = o (1.10c)

The successive overrelaxation method and symmetric overrelaxation

method are both slight modifications of the Gauss-Seidel method.

Two methods frequently used in common are the Secant method for
simultaneous nonlinear equations developed by Wolfe and Phillip[6].
and quadratically convergent Newton-like method based upon Gaussian
elimination deveIOped by Brown[7]. The two have been implemented in
code as ZSCNT and ZSYSTN respectively and collected into the INSL

subroutine liberaryl8].

In Wolfe's method(ZSCNT) at each step of the iterative process

there are n+1 trial solutions x1 .x’ , .... x(n+1) . [ultipligrg ..

J I

j=1.2. ... .n+1 are determined by solving the linear system

.j = 1 (1.11)
jg!
mi:

2 ajfi(xi) = o 1 = 1.2.....n (1.12)
j=1

Then the new trial solution is defined by

n+1
X =21j1j (1.13)
jg:

The Brown's method is based on Gaussian elimination in such a way
that the most recsent information is always used at each step of
algorithm. The modification suggests linearization of the components
sequentially. using each linear equation to eliminate a single compo-
nent of the solution from the remaining non-linear equations. as in
Gauss elimination. The system eventually reduces to a single
non-linear equation in a single unknown to which one step of the New-
ton iteration is then applied. The new values of all eliminated
components are then obtained in the reverse order by back substitue

tion.
1.3 Bond Graphs and Algebraic Loops

The existence of algebraic leaps in the equations of a physical
system may not be detected until the sorting or reducing process
starts in most traditional simulation approaches. But their existence
can be verified even before equation formulation when the bond graph

approach is employed. Bond graphs. which are based on energy storage

and power flow. allow system analysts and engineers to construct
models of electrical. magnetic. mechanical. hydraulic. pneumatic.
thermal and other systems using only a rather small set of ideal ele-

mentsl9].

The functional nature of the parts of a bond graph model can be
classified into the source field. energy storage field and dissipation
field. while the manner in which the parts interact can be represented
by the junction structure. All of this is done in a graphical format.
A bond graph model and its key vectors may be represented schematical-
ly as shown in Figure 1-1. The key vectors are labeled on each arrow:
U is the input to the junction structure from the source field; V is
the output from the junction structure to the source field; 2 is the
co-enersy variables vector; 1 is the time derivative of the energy
variables vector; Do is the input to the junction structure from the
dissipation field (generally a mixture of efforts and flows) and Di 1.
the input to the dissipation field from the junction structure (gener-
ally a mixture of efforts and flows). From the fiqure we see that X

is found from the junction structure as follows

Me

' “£(Z: Do- 0) (1.14.)

”e

or - 81,2 + 8,,00 + 81,0 (1.14b)

if TF and GI elements are constant.
For the storage and dissipation field vectors. we have

2 = of(X) (1.15)

 

 

 

 

 

 

 

 

 

 

SOURCE FIELD
( SE. SF )
U V
.—_—1| I ...—
l_——_I !—————
(C. I) I JUNCTION STRUCTURE Di ( R )
( 0. 1. IF. GY )
F 2 Do L
L_______
STORAGE DISSIPATION
FIELD FIELD
Figure 1-1 Bond graph structure and key variables

 

 

 

and U = U(t) (1.17)

for the source field.

From the above schematic diagram . D1 , the input to the dissipa-

tion field from the junction structure is given by

Di = GL(Za Do. U) (1.18‘)
or Di - 8,12 + 3,,00 + s,.u (1.18b)

if the TF and CY elements are constant.

Combining Eqs.(1.16) and (1.18b) we get

Di = 8312 + 813%(01) + sap" (1.19)

Depending upon the interaction between 8,, and 0L it may or may
not be possible to solve for Di from Eq.(1.l9) and then find Do. Tb
illustrate a case let us first consider a physical device shown in
Figure 1-2(a). In Figure 1-2(b). the corresponding bond-graph model
has been built. The I element represents the inertial effect and the
compliance effect is indicated by a C element in the mechanical sys-
tem. The R elements represent energy dissipative effects. The SF
element indicates an imposed velocity on the left (massless) plate as

an input.

Associated with every bond are two power variables --- effort and

flow. which are force and velocity in this mechanical device. respec-

 

 

 

 

 

 

 

 

 

L, I.
( a )
I I
1 v1 1 V1
3 4 ...—LS-Of—ES-é-
say—vo—=-o——=-nb SF Vo Kb
5 5
1 1
7 6\ / v6
6
C Rh C R;
( b ) ( c )

Figure 1-2 erield in a physical system
(a) Physical system
(b) Bond graph model with acausal bonds
(c) Bond graph model with causality assigned

11

tively. In this example there are six bonds; hence. there are six
force variables and six velocity variables. In addition there are two
(energy) state variables. p and 6. representing the momentum and the
spring deflection. respectively. There are 14 equations imposed by
the bond-graph structure through junction constraints and the field
constitutive relationships. It is desired that the state-space equa-

tions be obtained in an explicit form as follows

p1 = “1(p10 63. V.) (1.20.)
83 = “3(p10 5:: V.) (1020b)

where superdot denotes a time derivative.

Causality can be assigned to the bond-graph of Figure l-2(b)
according to the general rules[9]. After finishing the first step
(assigning required causality to the source SF) and the second step
(assigning the integral causality to the storage elements C and I). we
find that the causality does not fully extend through the graph. Some

acausal bonds (bonds 4.5.and 6) will be left (Figure 1-2(b)).

At this stage. we realize that an R-field exists in this system.
This implies that there will be an algebraic 100p in the system equa-

tions.

Suppose we continue the causality assignment by imposing an
arbitrary causal orientation on one of the two R elements. say. R‘,

Then we extend the causal implication through the graph using the con-

12

straint elements (0's and 1's). Now the causality assignment has been
completed(Figure l-2(c)). The state vector 1 and input vector U are

identified as follows

p; U = I V. l (1.21)

6

If '0 define F. and V‘ as auxiliary variables. then the system

equations are

i, s F. (1.22a)

8, = v. (1.22s)
and the constitutive equations are

Fe ‘ 34(va) ‘ Ia‘vo ' P1/I1 ’ V‘) (1.22!)

V. ' s.(F.) ' s.(F. - k5. ) (1.22s)

Assume that both R. and R. are linear.that is.

F4 ‘ R‘V4 (1.23.)

v, . a;‘p. (1.23s)

After some manipulations to eliminate the auxiliary variables F.

and V., an explicit state-space equation set can be developed; namely.

fl. = -[R.R./(R.+R.)m]p1 + [R.k/(R,+R.)]6a + [R.R./(R.+R{‘})]V.
(1.24a)

13
Pa = “[34/(R.+R.)m]p1 - [k/(R.+R.)]6, + [R./(R.+R.)]V.
(1.24b)

Now. 3uPP03° thIt the R. and R, are non-linear. Then we may have
difficulty solving the auxiliary equations to get an explicit state
form. In general explicit analytic solutions of non-linear coupled

equations are difficult. if not impossible. to achieve.

From the development above. we see that the process of causality
assignment is an aid in the process of identifying the algebraic loops
in dynamic systems. Furthermore. algebraic loOps in the mathematical
sense are physically related to the existence of dissipation fields.
Reading the partial causally-assigned bond-graph. we easily can iden-
tify the R-fields from other (dynamic) fields. This work can be done

by a digital computer automatically.

1.4 An Approach Stressing Parallelism

An approach to the simulation of this kind of systems containing
non-linear algebraic loops is the bond graph augmentation methodIlO].
The basic philosophy is to convert an algebraic loop subsystem into a
dynamic subsystem that conserves the intrinsic static characteristics
at its steady state and to employ a two-time-scale integration techni-

que.

This approach divides a large coupled nonlinear system into a

14

group of static subsystems and a group of dynamic subsystems. Every
subsystem in the static group is independent of the others. as is
every subsystem in the dynamic group. However the members of one
group will communicate with the members of the other group. At each
global time step. the static characteristics of every augmented sub-
system at its steady state can be computed independently by the same
algorithm. Therefore they may be done concurrently or in parallel.
The parallelism also can be applied to the computation of the dynamic
performances of the dynamic group at each global time step. The
two-time-scale integration technique places special emphasis on the
parallelism. which will become more practical with advent of the par-

allel computers.

The method includes mainly the following steps:

1) lake a bond-graph model for the non~linear system;

2) Identify and isolate the algebraic leaps before equation
formulation;

3) Dynamically augment the isolated static subsystems;

4) Select the proper parameters for the dynamicizers;

5) Use a twa-time-scale integration scheme to get simulation results.

In the following chapters we will discuss the details of this

approach.

15

2 BOND GRAPH AUGNENTATION lETHOD

2.1 Identification and Partitioning of erields

A bond graph gives a tep010gical picture of how the power (or
energy) flows in a particular physical system. Since every band con-
tains the power variables (effort and flow) and the causality shows
input/output relationship between the element pair. we can identify
R-fields readily. Let us consider an example. The bend-graph model
is shown in Figure 2-1(c). which represents a mechanical system (Fig-

ure 2-1(b)) or an electrical system (Figure 2-l(a)).

After assigning causality to the source element SF and energy
elements C and I. an R-field has been revealed. In this example . to
the right of bond 12 is a causally complete segment. while to the left
is a causally incomplete segment (Figure 2-2). The two segments share
bond 12. If we break the bond-graph into two subgraphs at band 12 .
we see that at any instant the output of the dynamic subsystem
represented by the right part is just the input to the static subsys-
tem indicated by the left part. and vice versa. The subsystems
interface at band 12. For example . at time t=t., en 1; obtained
from the dynamic field and input as an effort source to the static
field. If the output of the static field can be calculated in some
way which will be discussed in the later section. then the output will
act as a flow source at time t=t+AT to the dynamic field. by integrat-

ing one time step. AT. the dynamic field produce a new output to the

16

 

 

I(t)¢

 

 

L

( a )
l-—_V(t)
II I'M jl [’44
__JVZJ_JF77I c _____

R R R R R

( b )
(8) (9) 9 (10) (11) (12) (13) (14)1 (7)

SF'LOdILoLlr—z‘s—H 13 —H1—-11

R R R R
(1) (2) (3) (4) (5) (6)

( c )

Figure 2-1 Physical systems with a coupled R-field
(a) Electrical Circuit
(b) Mechanical system
(a) Bond graph model

17

1. . _
all all: all: mm.“ ..mwiml “14... (file (loll:

cu

asou-huaa- memouauuuam o:H NIN eunnmm

n_omu ems-aha uuoau nan-am

w n v m N
J 11!J

a

m

a

18

static field again. In this manner. the simulation of the global sys-
tem may be completed by interaction of its two subsystems as long as

the static field can be manipulated easily and efficiently.

In the general sense a global bond-graph can be partitioned into
causally complete and incomplete segments. The causally complete seg-
ments generally are comprised of energy elements (C and I).
dissipative elements (R). input elements (SE and SF) and the junction
structure (0.1. TF.GY). while the causally incomplete segments contain
the dissipation fields. input fields and associated junction struc-

ture.

In Figure 2-3 the concept of a general partitioned bend-graph is
illustrated. The interaction between the i-th dynamic subsystem and
the j-th static subsystem can be defined in vector notation. Each
subsystem may be viewed as an independent system with both input and

output vectors ascribed to it.

2.2 Dynamic Augmentation of Static Subsystems

The Basis of Augmentation

The basic idea in solving the algebraic loops is to convert a

coupled static field into a dynamic subsystem by introducing some

dynamic effects. to find the steady state output under a set of con-

stant inputs. and then to interact with other subsystems in the global

19

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Sn

 

 

 

fi1
e,:(f.i) 'D
-» j
fij (eij)
f
- Dn

 

 

 

 

 

 

Figure 2-3 The interactions among the
partitioned subsystems
Si : the i-th static subsystem
D5 : the j-th dynamic subsystem

20

system.

The dynamic augmentation can be done by introducing I- and C-
elements in such a way that the causality of the static field is fully

completed and the boundary conditions are not changed.

Consider a static field and a preposed dynamically augmented sub-
system associated with the static field (Figure 2-4). From the
definitions of the junction elements the following statement can be

made:

E efforts = O l-junctien (2.1a)

2 flows = 0 O-junctien (2.1b)
for the static field.

E efforts = p. l-junction (2.2a)

E flows = q, O-junction (2.2b)

and for the augmented system. Note that p1 and q, measure the energy

in I and C. respectively.

The condition of steady state requires that

Consequently at steady state the state equations will be identi-

cal to the junction constraint equations in the static field; i.e.

21

l I

R R

(a)

C [q] I} [p]

1 I

SFI—- 0—H1|-——-sr-:
R R
(b)

Figure 2-4 A static field and its augmentation
(a) A partitioned static field

(b) The augmented subsystem

22

E efforts = 0

2 flows = O

From this result it is evident that the algebraic character of
the static subsystem is preserved at the steady state of a preperly
augmented dynamic subsystem. At steady state the output c1 and f, of
the subsystem which are treated as inputs to the adjacent subsystems

can be easily determined by

91 = (IIC)q (2.3a)

f; '(1/1)p (2.3b)

To understand the nature of dynamic augmentation. let us inter-

pret the band graph as a physical device as shown in Figure 2-5.

Consider Figure 2-5(a). where A is a massless element being acted
upon by a force F. a velocity source V through a damper. and by anoth-
er resistive force from the ground. This really is a static system.
If a mass effect I and an elastic effect C are added into the system .
the sum of the forces at any instant is equal to the change in momen-
tum of the mass. and the sum of the velocities at any instant is equal
to the change in the length of the spring. At steady state. the rates
of change of momentum and length are zero; hence the dynamic system

represents the unaugmented static subsystem at equilibrim.

Only adding I elements to 1-junctions and C elements to

 

Figure 2-5

23

Schematic diagrams
(a) Static system

(b) Modified system

 

 

(b)

N

24

O-junctions facilitates the completion of the causality assignment and
does not change the original causality orientation at the boundary

bonds. So. the principle should be

(1) Add I elements to the l-junctions

(2) Add C elements to the O-junctions

The added I and C elements are assumed to be linear. conservative

energy storage elements.

Some Definitions

For convenience in the development later on . it is useful to

define some terminology which is used frequently in this chapter.

Definition: An erield or static field . is a collection of dis-
sipative elements (R). source elements (SE and SF) and junctions
(0.1.TF and GT). in which the causality assignment is not able to be

finished uniquely.

Definition: A dynamic field is a segment of the global
bond-graph in which there is at least one dynamic element and the cau-

sality can be assigned completely.

The dynamic field may be a collection of some basic elements

25

(C.I.R.SE.SF.O.1.TF.GY) or it may contain only one single dynamic ele-

ment (C or I).

Definition: An intgrttting‘hggg is the connecting bond between

an R-field and a dynamic field in the global bond-graph.

Definition: An acting junctton is a junction with one or more

interacting bends attached.

Definition: A gglti-ttangh tng juggtion is a junction which

possesses more than two internal bonds.

Definition: A tingle-R type jungtion is a junction to which only

one R element is attached .

Definition: A mglti-R typg jgngtien is a junction to which more

than one R element is attached.

Definition: An tmtginary tource is a constant effort or flow
source whose magnitude equals that of the output of the adjacent

field.

An example is shown in Figure 2-6. The global bond-graph model
has been partially assigned in causality and then partitioned into two

dynamic fields and one R-field. Bond 8 connects dynamic field 1 and

26

at :1 vi m1 1,

81 f DI
SE_‘-il"_~0—__- 1——-0——-'SF
l l l . SE—HIP—AR
a R a
l
I

Figure 2-6 Partitioning
(a) Global system

(b) Partitioned subsystems

27

the R-field. while bond 11 connects the R-field and dynamic field 2.
All the interactions among these subsystems occur on these two bands.
so they are the interacting bonds. Nodes (9). (11) and (15) are set-
ing junctions because interacting bonds 8 and 11 attach to them.
Nodes (6). (9) and (11) are single-R type junctions. while node (4) is
a multi-R type junction. In the partitioned subsystems. SE12. SF21.
SE11 and SFll are introduced to the subgraphs for graph completion and

input functions. They are defined as imaginary sources.

Order of the Augmented Subsystem

It is evident that every added energy storage element will intro-
duce a new state variable and a corresponding parameter. Hence it

will increase the dynamic order of the subsystem by one.

For a large R-field with a number of non-linear resistive ele-
ments. by increasing the number of added dynamic elements. the dynamic
order will be increased considerably. This increase may benefit us in
one way. in that it increases the flexibility in satisfying the
desired causality orientations of all R elements in the subsystem.
However. the computational efficiency may be decreased with the
increase of the dynamic order. Furthermore. the higher the dynamic
order. the more difficult parameter selection will be. because an
imprOper selection of parameters for the added elements may cause a
considerable delay in the field reaching its steady state .

Consequently augmentation greatly affects the computational efficiency

28

and the accuracy. Parameter selection will be discussed in detail in

the next chapter.

Augmentation could be done according to at least two different
considerations. One of them is to make the causality orientation
satisfy the nature of every particular R element. Reconsider the sub-
system containing an algebraic loop represented by the bond-graph in
Figure 2-7(a). By adopting the principle mentioned before. we have
three possibilities to augment the subsystem by (1) adding only one I
element. (2) adding only one C element. or (3) adding one I and one C
element to each 1- and O-junction. respectively. If the constitutive

functions of the resistive elements in the system are

01 3 31(f1)

e, - a,(r,)

choice (1) satisfies the required nature of these R elements. The
causality would be assigned as shown in Figure 2-7(b). Other choices

are shown in parts (c) and (d).

In this case there are two single-R type junctions; one is the
O-junction and the other one is l-junction. The maximum number of
useful added dynamic elements is two; that is. it equals the sum of
the number of O-junctions and l-junctions. However. in a general case
with more single-R type or multi-R type junctions. it could be visual-
ized that a number of possibilities to introduce the storage elements

exist. The number of useful ways depends upon the type of resistive

29

set—- 0 -——- 1 |-—-sr~:

SF}—. 0 —._-' 1,—555

C I
i a
1’}, C I

e B R (f
e I R (f

choice 1

f = G (e

f = G (e

choice 2

f = G (e

e = R (f

choice 3

e = R (f

f = G (e

unrealizable

Figure 2-7 Causality orientations

(a)

(b)

(c)

(d)

(e)

30

elements in the algebraic loop. If some of them are non-linear and
their preferred causality is required. then the number of ways in

which the C and I elements are selected is reduced.

Suppose that the causality orientations of all R elements in a
nefield are not of concern. That is. the inverse expressions of dis-
sipation functions are easy to get. Then the causality orientations
could be any combination. But the maximum order of augmentation would

be sum of the numbers of the O- and l-junctiens.

Carefully studying the example above(Figure 2-7) reveals the fact
that only three causal arrangements are available through selective
dynamic augmentation. However. the dissipative elements may exhibit
four different causal arrangements. For the fourth situation.
part(e). the method of dynamic augmentation is inadequate to force the

preferred causal orientation on both R elements.

Generally. for translational and rotational dissipative elements
and for most fluid resistive effects we have the case in which a dis-
sipation is the function of flow. But some cases. such as
semiconductors in electronics and fluid leakage in hydraulics and
pneumatics. have the inverse dissipation function. Fortunately. many
non-linear dissipation functions are not too complicated and their
inverse forms may be found without much effort. But important excep-
tions exist. For example. in mechanical systems the dry friction is

not allowed to be a function of force.

31

Augmentation Procedure

Assume that every dissipation function can be expressed either as
e = R(f) or f = G(e). Then it will be possible for a computer to aug-
ment the subsystem automatically without regard to the preferred

causality orientations. The procedure is designed as follows.

First we note that every subsystem needs to compute the interact-
ing bond variables as output to the adjacent dynamic subsystems.
These bend variables as outputs always are the common forces or flows
associated with the O-junctions or l-junctions. Therefore. if we add
energy elements to the acting junctions. we not only facilitate the
causality assignment but also obtain the output vector from the
co-energy variables without invoking more complicated algebraic output
equations. So adding the I or C elements to the acting junctions

possesses the highest priority.

For example. the subgraph in Figure 2-8(a) represents a static
field containing a real source SE1 and an imaginary source SE18. Bond
18 is an interacting bond and node (18) is an acting junction. At a
certain 810511 tifl° instant. 9;. is the input from an adjacent dynamic
field which holds constant during the integration in the local time
scale.'hen the augmented subsystem reaches its local steady state. f1.
will be the output to the dynamic field. This output is just the com-
mon flow on the acting junction. If we add an I element to the acting

junction (18) (Figure 2-8(b)). then the output can be calculated from

SE —-11—-o—-1—‘1-—15-1-15-o—11-1l——-sn

I /I\ I\ I ”I “I ”I

R R R
(a)

I
se—-l1-—-o—-\ 1\—o—1—o —-l1l-—-SE
I /I .I .I I I I
R R R a R R R
(b)

c /_I [I
sc-—-t

PEN l l l

R

. f f /‘
ss_qi|.._..'ol_.{ _.‘1|.__..o_..1,_.se
R

l/lxi'fj . 1

(d)

Figure 2-8 An example
(a) A static susystem
(b) After step 1

(c) After step 2

(d) After step 4

33

f18 = fie ‘ Pas [‘19

Second. the bonds associated with the multi-R type junction can
not be expected to complete a causal assignment by extending from
neighbor junctions. We complete the causal assignment by adding a
dynamic element directly to this type of junction. In Figure 2-8(c).
a C element has been added to node (12) and an I element has been
added to node (13). In a computer program it should be done in des-

cending order until no multi-R type junctions exist.

Third. if there exist some multi-branch type junctions in the

static field. it is desireable to add dynamic elements to them.

Fourth. add the dynamic elements only to the remaining unassigned
single-R type 0- or 1-junctions. whichever is fewer in number. If
these two kinds of unassigned junctions have the same number. either

kind can be chosen for augmentation.

In this example. there are two single-R type O-junctions and one
single-R type I-junction unassigned. For the lower order augmentation
adding one I element is enough to complete the causality assingment

for this subsystem (Figure 2-8(d)).

After assigning and extending the required causalities for all
the sources and the preferred causalities for all the storage

elementes (C or I). the program PART will begin to identify the dynam-

34

ic fields and the static fields. then partition them into subsystems
by using the imaginary sources to preserve the boundary conditions.
All the information will be stored in matrix form where every column
or row contains a dynamic or static subsystem. The program ADGNEN
finds the most appropriate 0- or 1-junction to add the corresponding C
or I element according to the procedure discussed before. They will
be incorporated into the ENPORT-G program. The calling tree is shown

in Appendix A and the program listings are given in Appendix B.

2.3 Two-time-scale Integration

The solution for the steady-state vector X in the case of linear
dissipative fields can be achieved by simple linear algebra. provided

the local [A] matrix is nonsingular.

x = -[A]"[B]U

However. in the case of non-linear dissipative fields we can not
use this approach to get a solution. Although we can employ a linear-
izatien method to get a constant [R] matrix at a certain time instant.
then use the formula above to calculate steady state. it will reduce
the accuracy of computation. For example. we can calculate X at time

t through

xt = ’[Alt-Atwwt

35

The [A] matrix is not the present value. If we can use
non-linear integration to get the steady-state in the local time
scale. the computation result will be more accurate and so will the

output vector.

The augmented subsystems need to reach steady state by integra-
tion under a set of constant inputs which are the outputs from
adjacent dynamic subsystems at global time T . It seems that when the
dynamically augmented subsystem is under the process of integration in
local time . t . the global time is 'frozen'. After the steady state
has been reached. the output vector of augmented subsystem will be
sent to the adjacent dynamic subsystems and the global time will warm
up and make a step (AT) forward for the global integration. A new
output vector from the one-step integration refreshes the input vector
of the augmented subsystems. These subsystems will integrate again in
their local time scales to get a new steady-state . This process is
repeated for the total global time interval and the desired simulation
results. Since there are two time scales in the whole simulation pro-
cess . one for the local subsystems and the other one for the global
system. we call this scheme two-time-scale integration. A schematic

diagram is shown in Figure 2-9.

36

 

I Global systexT]

T

 

 

I

,LI Input e1 (or fij’
I into the subsystem

t I

Locally integrate to find
steady state for the
dynamically augmented subsystem

 

 

 

 

 

 

New input eij(or fij)
input into the subsystem
at T+dT

 

 

 

Output fji(or °ji)
to the global system

_ I

Global systemintegrated
over time step dT

 

 

 

 

TWdT

 

 

Figure 2-9 Schematic diagram of two-time-scale integration

37

3 SOME COMPUTATIONAL CONSIDERATIONS
3.1 Effect of Eigenvalues on Computation Efficiency

Our ultimate goal is to simulate the dynamic response of
non-linear systems. but the following discussion on linear systems may

benefit this.

For a dynamically augmented subsystem. the system representation

is readily arranged into an explicit state-space form

i = [Alx + [BIU (3.1)

Structurally. the [A] matrix may be resolved further into the
following form since no dependent energy variables exist (i.e.. all

integral causality):

[A] = -[ a ][ u” ] (3.2a)
for the augmented subsystems by adding only I elements.
or [A] . -[ a ][ x ] ' (3.2a)

for the augmented subsystems by adding only C elements.

The [R] and [6] matrices are derived from the bond-graph topology
and dissipative elements. while the [IT’] and [R] matrices consist of

the free parameters introduced through the dynamic augmentation.

38

The [IT‘] or [R] matrix corresponding to the energy-storage
fields is diagonal and positive definite for integral causality. all
non-zero terms in the [N"] or [R] matrix are positive. For the
Refields to dissipate energy for any possible pert condition. the
resistance parameters should be positive . This is true because only
the 1-port resistors can dissipate power and the junction elements
conserve power. From the preceding statements. it follows that the

[R] matrix will be positive definite.

The question that arises is ”How should the free parameters be
selected to provide computationally efficient convergence to the

steady state?"

As mentioned before. during the local time integration the input
vector is a constant vector. We want the augmented state variables of

each subfield to follow a path like below to reach their steady state

x = $1: eAt k. {5" BO a: (3.3)

where A is an n n matrix with real negative eigenvalues and eAt is the
fundamental matrix of the system. B is a n‘m matrix and U is the input

vector with dimension m.

Let us call this form "diffusion"-type convergence. Graphically.
the "diffusion" convergence has the exponential feature as in Figure
3-1(a) and 3-1(b).The n‘n matrix [A] has n eigenvalues. Suppose that

the augmented subsystem is a ”diffusion"-type. that is all the eigen-

39

0.673%

 

 

N

.. /

 

 

 

 

 

 

 

0.632
X
i ’_'-—
lls 2/s 3/s 4/s 5/s t
(a)
X
xi
Xi ‘_4_
O lls 2/s 3Is 4Is Sls 4=;—t
(b)

09
Figure 3-1 Diffusion" convergence

40

values are negative real numbers. They are s1, .3, ..., 'n' To
achieve reasonable accuracy the integration duration should be set to
at least five times ( 1/I‘Imin) and the time-step should be set less
than one third of 1""max (Figure 3-1). From these considerations.
we hepe all eigenvalues are clustered. This means that all the eigen-
values lie within a amall region in the s-plane. so that the number of

integration steps can be minimized .

Mathematically. we can state the objective as trying to choose
kii such that the normalized ”spread” of the eigenvalues. p. be minim-

ized. namely

9 g ( [‘Imax ' I‘Imin ) I ( I‘Imax + [‘Imin ) (3'4)

where ‘min and ‘max are the minimum and maximum among the all

eigenvalues.respectively. The value of the "spread” varies between

zero and one .

The eigenvalues of the [A] matrix are determined by the charac-
teristics of the [R] matrix. which is intrinsic to the field . and the
[R] matrix. in which every diagonal term is the parameter of an added
dynamic element. Although the minimum of the ”spread" p is determined
by the [R] matrix itself. getting to the minimum depends upon how we

select the parameters of the added elements.

41

3.2 Parameter Selection in Linear Dissipative Fields

Now let us discuss the issue of how to select the parameters of
the added elements starting from the simplest first-order case to get

same feeling about the general rule for the n-th order case.

First-order Augmented Subsystems

The first-order augmentation has one free parameter to be select-
ed. It is the simplest case in Refield augmentation. Consider the
R-field shown in Figure 3-2(a) and (b). The state equations for the

augmented erields (a) and (b) are given by Eq.(3.5a) and Eq.(3.5b).

respectively:
i‘ = -(Rz + R‘)p‘/m‘ + szl - F. (3a5‘)
i. - -(n. + n.)p,/a, + p, - s, (3.5a)

where the [A] matrix is

[A] = ’[R3 + R4](l/m) (3.6)

The eigenvalues of both equations have the same form

s = -(R3 + R‘)/m (3.7)

If we set s to be an arbitrary constant. n . for example. s = -1.

the parameter m may be determined by

42

I
SFr—-o—-- 1|-—-—sa sat—l—o—i-{ié‘ss
I I at I
a a a a

( a )

ss—q1}——-ss SE—l—lll-L-SE
/'--\ 7 ‘
H R R R
(b)

Figure 3-2 The first-order augmentation

(a) Single-R type
(b) Multi-R type

V1 [p10] V1 [Pu]

1 3 s 7 9
SF }———-- o'--‘i 1’F-—“ 0 """i 1 ["""SE
2 4 6 8
a a a R
p s
x ‘ 10 U . 1
e
P11 9

Figure 3-3 The second-order augmentation

43

n = R. + R. (3.8)

Using the eigenvalue s c -1. we can choose the time interval and

the duration for the local time integration

At = (1/3)‘(1/s ) 8 1/3 (3.9a)

and tf = 5‘(l/s ) = 5 (3.9b)

which requires 15 steps for "reasonable" accuracy.

Second-order Augmented Subsystems

A typical second-order dynamic augmented subsystem is shown in

Figure 3-3. The state equations in matrix form are

i.. -(n. +3. +2.) a. /n,. o p.. a, 0 v1 (3.10)
I +
$1, a. —(a. +a. ) o 1/nu p,, o -1 F,
where the [A] matrix is
[A] ' -(Ra +34 +Rsumac Rs/nai
Rs/nie ’(Rs +Rs’llaa

At this point it is the appropriate time to apply the Gershgorin

circle theorem to study the eigenvalue distribution.

The Gershgorin circle theorem[ll] states that:

44

Let [A] be an n‘n matrix. and let Ci, i a 1, 2, ..., n be the
n
discs with centers a-. and radii Ri= 2 ‘ik . Let D denote the union

11
k=5kfii
of the discs Ci' Then all the eigenvalues of [A] lie within D.

If we inspect the [R] matrix. which is generated from the
second-order linear system. we may find the off-diagonal terms are
less than the corresponding diagonal terms. This implies that the
locations of the circle centers are more important than the radii in
estimating the distribution of eigenvalues for this type matrix.
Generally we expect that if all the circle centers are located at the
same point in the s-plane and all the radii are as similar as possi-
ble. then the distribution region would be reduced. Therefore. the

"spread” of eigenvalues may be the minimum.

For simplicity. let all centers be at -1.0. The conditions are

Rs +34 +3s ‘ Inis (3.11m)

Rs +Rs ‘ '11 (3.11b)
then the [A] matrix becomes
[A] '3 '1 Rg/ (Rg +ng)

R./ (a, +3. +a‘) -1

The characteristic polynomial p(s) would be

45

p(s) = det(sI - A)

= s‘ + 2s + 1 - R:/[(R,+R.+R.)(Rs+la)l ‘ 0 (3°12)

The roots of p(s) are

31,, = -1 3; [xi/(aga‘an(a,+a,)]‘/’ (3.13)

Because the radical term is always greater than zero and less
than one. two negative real eigenvalues are produced. Therefore

”diffusion" dynamics could be realized.

We also can analytically verify the correctness of applying the
Gershgorin circle theorem using a general second-order augmented case

in which the [A] matrix has the following form

[A] T T 311 81a ka 0

“as 8a: 0 ka

where the [6] matrix is symmetric. positive definite and the [R]
matrix is positive diagonal because of the topology characteristics of

R-field. Another preperty of matrix [6] is that the diagonal terms

811 are dominant.

The characteristic polynomial p(s) would be

46

p(s) = det(sI - A )

' ‘1 + ‘(Uiaki + Beaks) + (RiakaBaaks ’ Riskaks) (3°14)

The roots of p(s) are
81,. = -(s..k. + leaks)/2 :.l(s..ki + saak.)’
- 4(s,,k.g,,k, - ;:,r.k,)]"’/2
(3.15s)
°’ ‘1.2 ‘ ”(811‘s + s..k.)/2 i.l(s..k. - s..k.)’ - 4si,k.k.]‘/’/2

(3.15b)

Eq.(3.15) shows that the roots will be real and negative. We

simplify the notation of s1 : as follows

.1,, = -a/2 1 (3’ — 4C)‘/’/2 (3.16)
where

B = a..k. + s..k. (3.17)

C a 8111.33.11 ' SIakaka (3°18)

Then the "spread" is

p = (1 - 4CIB’)‘/’ (3.19)

Since p > 0. we can minimize p’ with respect to k1 and k,. namely

de’)/dk. = —4(a.....k.-::.k.-2c:../B)B’
(3.20s)
d(p‘)/dk3 = -4(g,,g,,k,~g§,r,-2c3,,/B)a“

(3.20b)

47

These equations yield

2Cg11 - BC/k; = 0 (3.21s)
and 2Cg33 - BC/ka = 0 (3.21b)
Finally we have

kaiSii ‘ k283: (3.22)

To ensure a minimum we also can show that d’(p')/dk: ) O and

d’(p’)/dk: > 0 . So if we choose k, = 1/g11. then k, would be 1");

from the above relation. i.e.. one of the optimal selection will be

ks ‘ 1,811 (3.23a)

ks ‘ 1,83: (3.23b)

This coincides with the result obtained by using the Gershgorin
circle theorem. Let k1 = l/g1; and k, I 1,833. then "spread” p will

be

p = s..(1/g,,g,,)"’ (3.24)

Since 3:3 < 8113:3' ‘0 P: T “Ia/81132: ( 13 hence P < 1 °

The above observation has been illustrated by an arbitrary numer-

ical example in TABLE 1. The [A] matrix is

48

 

 

-----l"-

 

ao=~¢>commo

 

 

 

 

 

 

 

 

 

 

 

 

 

hamm.o «aaa.o cfiwc.o mo vacuum
wee.ot was.cu vee.cu Na
. mo=~a>commm
eoo.~u moo on- eec.~I am
cNN "Amu+e=vu um uAmz+wmvcnxu c—~ nmx+wm an: aumoaogo momma
o—N umm+vm+~= can uwm+v¢+~m saw uwm+vm+uz Na; no amouosmumm
oc— m:
cm um «amuse—o a ma
can we «Homoseuum
can u“
m mm<u N mw<u u mm<u
ZO~P<hzmzc=< «mono czcumw < no mmag<>zucum «In mac<b

49

[A] - -2io/au 10/.n

10"10 -110/I11

Three different parameter sets of In and mu have been used. Their
eigenvalue distribution regions are shown in Figure 3-4. Obviously

the parameter pair in case 1 e the optimal selection among the choices

shown.

General n th-order Augmented Subsystems

Experience indicates that many physical systems containing alge-
braic loops may be adequately handled with only first or second order
augmentation. But it is still possible for higher order augmentations
to emerge in some complicated systems. A useful rule may be general-
ized from the discussion above on the first and the second-order cases
to the general nth-order case. In general. the augmentation by adding
only I or C elements will produce the resolvent [R] and [N"] or [G]

and [K] matrices as follows

. 1 . .
Isa 113 ...... rgn III: ..... 0
r3: r,, ...... . . 1/m,

[A] -—ta][u“] = - . . . . .
rn1 ...... rnn LO ..... llmn

 

 

 

 

-1.00
0.048

-1.00
0.048

-1.00
0.048

50

 

 

 

 

 

AIm
0.091
-1.00
-1.0 -005 O fie
(a)
111m
0.909 ] l
-10.00 ‘ t1
E3-10.0 II -1.0 0*452
(b)
II I...
0.045 ]
-o.500
A A fl.-
-1f6' $015 0 Re
(a)

Figure 3-4 The eigenvalue distribution
(a) Case 1
(b) Case 2

(c) Case 3

 

51

or FRI, gu ...... Ian Ik1 ..... 0 .
‘3; ‘33 000000 a o k;
[A] = -[G][R] = - . . . . .
.‘n1 aaaaaa ‘nnJ .0 as... kn“

 

 

 

 

The general [R] or [6] matrix has following properties as men-

tioned before:

1) The [R] or [G] matrix is symmetric and positive definite;

2) The diagonal terms rii or 3ii are dominant among the terms in

the same row.

The idea that the free parameters should make the diagonal terms
an and a,, in the two dimensional product matrix [A] be minus one so
that the ”spread” of eigenvalues is the minimum. may be extended to a
general n-th order augmentation. So the general rule for optimal

parameter selection may be stated as that:

The Optimal parameters should make all the diagonal terms ‘ii in
the product matrix [A] be an identical number. more precisely. m1 =

‘ii for adding only I elements or ki ' 1"ii for adding only C ele-

ments 0

52

The most common junction structure of Refields is so~called
"chain" type 0/1 junction structure (Figure 3-5) in which the junction
elements 0 and 1 are arranged in alternation. It may be worth paying
a little more attention to it. Obviously every added dynamic element
increases the number of eigenvalues of the [A] matrix by one. The
additional eigenvalue will be determined principally by the parameter
of the dynamic element itself and the values of the dissipative ele-

ments in the neighborhood.

The [R] matrix produced from a general n-th-order dynamic
augmented static field with 'chain' junction structure by adding I or

C elements possesses a banded form:

r
”(R:+R3+Rg) R. 0 O O O O
R, —(n,+3‘+R,) R, 0 0 0 0
-[R] g 0 e o o o a o
e e e Ri_1 -(Ri-1+ni+ni+z) 31+: 0
0 ° ° ° ' Rn-1 ’(Rn-1+Rnxl
L

 

An explicit rule of parameter selection for "chain" type junction

structure may be stated as

 

53

If“
ssh—.0 —.o_..11}_.0—- 1|——- 82

I III I

R1321-1 II21 ll21+1 Rn

Ii " R21-1 + R2i "’ R21+1
( a )

(lei

sat—*0 ——-l1I——0——-11I——- 1}—-ss

I Ill 1

R1321-1 R21.1 R21+1 Rn

C1 = IIRZI'I + l/Rzi 4’ 1lR2i+1

( b )

Figure 3-5 General n-th-order augmentation
(a) By adding I-elements
(b) By adding C-elements

54

1) for augmentation by adding only I elements. use
'1 ‘ ‘21-: +3.1 +3.... (3.25.)
2) for augmentation by adding only C elements. use
1/ki = 1/R,i_1 + 1/R,i + 1/R,i+1 (3.25b)
where “i --- the m parameter of the i-th introduced I element;
k- ---- the k parameter of the i-th introduced C element;

1

1 ---- linear resistance in the field.

An example for general case is shown in Figure 3-6. The [G] and

[R] can be derived from the topology of the bond graph model as below

 

 

'l(1/k.+1/a.+1/a,+1/a.) 1/n, 1/n. ‘I I1.. 0 0'
-[G][x] - 11R, -(1/a.+1/a.) 0 o 1.. 0
I Ill, 0 -(1/a,+1/a.) L0 0 k1,

 

 

Some numerical results are tabulated in Table 2. It shows that

the optimal parameter set is obtained by using the suggested rule.

It has been noticed that for those [R] matrices in which some of

thO rij terms ( i - j ) are greater than zero. the eigenvalue "spread”
obtained by the general rule is not the minimum. However. employing
the general rule still makes the parameter set result in a satisfac-
tary eigenvalue distribution in the s-plane. Although the
mathematical proof has not been given. the general rule still can be

applied to such [R] matrices. An example is given in Figure

3-7(notice that this case is not of the least augmentation). The [A]

55

”I

1_5_.‘ass
12T1sc 1 c a“
s. \/. .. f1.
—-7-Ii]—-0—-IH1l-—-0—-|ss
11 2 31 4
R1 R2 “3 “a
( b

Figure 3-6 An R-field: example 1

(a) Before augmentation

(b) After augmentation

56

 

.uoa acausauen nesuuao ash

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

IIII a

~eas.c m~sw.¢( wmoM.OI nvom.ul o.m c.“ n.~

emce.o cmN~.~I ~h-.~) mmvm.o1 c.n c.cH m.~

-~h.c upmo.cl ccoo.—I wuom.vl o.m c.n c.cn cu an cu cu ca c.o~
a cccné seen. u) cocoJu ccOn.OI c.m c.m n.~
...... -... . 2r...

newn.o wwNn.ct oco~.0I unhv.cl O.“ O.~ a.“ fi

am um um hug bum max on ma v: mu «a an
I IIIIIII VII-l. I -LIIIIIFIIII-[fi [I 11111 fl» IIIITIIIILrIII

cocoon no=~e>cowum uuouoeeuem comm «amuse—o a mo euoaosauea

.IIIIIIIIIIILrIIIIIuIIIIIIIIIIII I IIIIIIIIIIIIII

 

 

 

 

W'l:-'- . ----

zc~h<hzmzus< umamo akm=h < no mu=4<>zuunm «In man<H

 

57

11 5,
Sis—‘41 a5
6 7 s 9
sa—-|1———~0 ——-1—--0i-|sa
1 2] 3 a

1 R2 3 4
( a )

14
ss-il-l 1l-L- a5

lez /‘

se—-|1|——7- o——-|1|——.0 -—°—-|ss
11' 21’ .1 41’
a1 a2 a3 a4
(‘b )

Figure 3-7 An R-field: examlpe 2

a) Before augmentation

b) After augmentation

58

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

«oom.o upp.ct— vnm.91 uhn.ut o.ov o.cw c.on
I. mnan.o11 (mo~.~1_..lmmn.nu awn.cw. o.c~. c.om . c.9n an c“ on on on >H
a mopv.e www.ct cmh.o1 nuo.ua c.c~ c.ow c.on
v—ww.c can.HI ovo.ml cc.nl o." o.« o.«
~awm.¢ vna.ci no.~1 vc.nt o.~u c.- o.nu
cnmn.c oo~.«t «cm.o1 gov.) o.«n m.n— °.~v om o" o« N ca muu
a pu-.o ~va.e) mom.¢l enu.nt a.nn c.nn o.uu
o~an.c owwnnu oo.o—I ca.-| o.« o." O.“
chaa.o hnu.01 n-.ol uvu.~) c.ow o.oc~ $.95
ecoa.c who.~1 pwu.ou ac~.oa o.co c.oe o.oo~
ream.o «en.~l www.01 ma~.OI o.°o o.oo o.co— an as an an c" an
Mllmmhm.o van.ul pu~.oi av~.cl o.oo o.ow c.oh
mvom.o n.nwn) cc.o—I ~n.m~t c." o." c.~
noou.o mnh.—I hco.ci anv.on c.c~ c.m~ o.cm
«vow.c vmo.~1 con.01 www.ct o.o~ o.o~ 0.0m
ennn.c onh.~n con.on eon.ou o.ou o.o~ o.ov an o“ o“ c" on H
a wnpn.o mna.u1 Huc.cl can.91 o.ou o.o~ o.On
meno.c onm.n«) oco.cul Nv~.vvn G.“ o.« c.u
mm um «m was vaz nu: ma cu m: n: as .
oaau
vacuum nee—aamouSm uuouoaauam ooum uumoao~o a he aneuoaauam

 

 

 

 

 

 

 

 

 

awhmhmuam ~¢<mhmaz< z< kc mm=4<>zmonu

mun mqnmh

 

59

matrix is derived as below

 

 

L(Rs'I'Rs'I'Ra) -R, as I1/m1, 0 0-
[A] = '3, '(R,+R,) Rs 0 1/m14 0
. R2 R3 -(R1+R3)J .0 0 1/m15J

 

 

Some numerical computational results on the eigenvalue ”spread"
are tabulated in Table 3. It shows the parameter selection by using
the general rule is rather close to the possible best one; therefore

it is acceptable.

The general rule for a "chain" type structure has been applied to
a number of arbitrary [R] matrices with different dimensions and
numerical conditions. Figure 3-8 represents a third-order augmented
R-field by adding three I elements which introduce three free parame-
ters: '1: m3 and m.. Three different cases have been
investigated(Tab1e 4). For the first case the parameter values of all
R elements are equal or close each other; we call this the normal
case. For the second case the parameter values of R elements are dif-
ferent in such a way that the R elements attached to the junctions to
which the dynamicizers are added have relatively large parameter
values compared with their neighbors. The third case is just the

opposite of case 2.

In every case. if we choose parameter values for m1, m, and m,

according to the general rule obtained before. then the "spread” p

60

11 12 13
SF}——-o—-|1}——-0—-11}—-o——-|'1I——-sa
R1 R2 R3 R4 33 36

Figure 3-8 A third-order augmentation

61

 

 

 

 

 

 

 

 

 

 

 

 

 

 

n: + v + n: a an

§+~z+§u~= Inna
«nom.o vco.ul ~n~.ct non.ui c.0a c.nm c.nn
huha.o nmc.o1 nvn.nt www.ct o.nv o.cn~ o.nn

oauo>m
«nee.o «no.0: «ov.~) oc«.on ounv c.nm c.0h m cm n on m c« <
~o«u.c cn~.o0 onm.~1 coc.~n o.nm c.mm c.nn
naav.o «ca.on vmv.oc vcu.nl c.cnn c.nr o.omn
coon.o com «I wnv.ot «mo.ut c.nm c an" a can om n on on ecu an o—a-uo>au
hnov.° nac.ut uwv.ol mma.°I c.nw o.nh c.om~
nvn~.o nu~.u) nah.o1 coo.n1 o.nw c.mh o.onn
nn—m.e www.ct vnm.cl nan.no c.cm o.on c.on
«can o man «I «on an coo «1 c on c cm 0 cm on as an on an cm macaw
nm¢o.° noh.o1 anv.nl www.cl c.c~ o.cm o.om
hmnn.c sun.n) mpv.cl ooc.ul o.o~ o.cm 0.9m
mm am um um um an m: n: v: n: N: a: aeouummeou

maoumm coma->monmm .ao—o momma uo .a-uam numoao—o mo uuouoaauam

 

 

 

 

 

NMPHO=MHm z~<=u =hmh awhmmmnam < he mm=4<>2uc~m min mqn<h

62

will be the smallest. Furthermore. among the three cases the smallest
one is achieved from case 2. Case 2 is called one with favorable con-
dition. Generally we hope the R-fields are in such favorable
condition which will make the local integration be more efficient.
The apposite case shows a much worse result. so it is called the case

with adverse condition.

3.3 Parameter Selection in Nonlinear R-fields

Linearization of R Elements

In a nonlinear dissipative problem the "parameters" of nonlinear
R elements vary. The eigenvalues of the associated [A] matrix change
during simulation. Consequently the parameters of the added dynamic
elements may need adjustment during simulation. Otherwise they may
cause integration in the local time scale to be unstable or lead to a

considerable error.

One technique for applying previous results to nonlinear system
problems is linearization. At selected times and states. the equi-
valent R-field parameters can be found by linearizing the dissipative

element characteristics. That is. we find the local tangent

R = de/df (3.27)

where e = s,(f). (3.28)

63

Alternatively a chordal approximation based on the global e and f

can be used. namely.

R = e/f (3.29)

For example. a nonlinear dissipasive R element possesses the fol-

lowing constitutive function

e = 0.1511-3 (3.30)

By using the chordal definition for the parameter of a linear R

element. the equivalent parameter of R would be calculated by

ac - e / f - 0.1s£‘-’/£ = 0.151°-‘ (3.31)
whereas the tangent approximation yields

R. = de/df = 0.1sf'-‘ . (3.32)

Reference R-Natrix for Integration Control

Partitioned static and dynamic fields are shown in Figure 2-2.
One of the possible augmented subsystem is preposed in Figure 3-9.
with the state vector. 1. the input vector. U. and the output vector.

V. The constitutive dissipation functions are

64

11 zl 3)
R R R
( a )
4;g;>;l [p10]
5
SFf—b 0 __-41'__- 0.___-'
1 {1' 3
R R R
x = U a: 5
p11 °9

Figure 3-9 A nonlinear R-field

(a) Before augmentation

(b) After augmentation

4l

R

1 |-—- sa
84-[

R

65

e. = 2.osraN(r.)£:°’ (3.33.)
e, = SIGN(1.)£: (3.33b)
e, = arana(r,) (3.33c)
e‘ = 2.01. (3.330)

The equivalent parameters of the R elements linearized by the

chordal approximation would be

R1 = 2.0310N(£.)£:-’ (3.34.)
R. = SIGNIfalf. (3.340)
a, = 3TANH(f,)/f, (3.34c)
a. = 2.0 (3.340)

The state equations for the nonlinear subsystem can be derived
from the bond graph model as usual. Remember that the added I ele-

ments are linear with constant parameters m1. and m11.

£1. - 2.0SIGN(r.-p,./n,.)(1,-p../n..)‘°' - SIGNIp../a..)(p../a..)’

- 3 T'ANH(p"/m1° - pat/m11) (3.35a)

3 TA”R(Pae/‘ao ' DIR/n11) ' 2°0p11/m11 ’ 99

p11

(3.35b)

If we arrange the above equations as follows

66

*o' (p,./n,.)’

(f - lm ) - ( /m )
(fl-plo’mxo) s Pie so (p‘./m;.) Pie in

2.0(fs’pae/mao)

 

 

p10

 

3TANH(p1./mgg - pit/n13)

(Pie/”as ' pan/M11) P‘./m1° - p1,/m11)

(3.36s)

3TANIi(p lm - p Im )
= 1. I. :: 1:(Pae/“io ' Pas/'11) ' Z'OPIII'II - o’
(pgg/mgo - Pia/mas)

 

p11

(3.36b)
Further manipulating yields
Pie 3 ' [2-0(fs ’ FRO/‘10)... + Pie/lie + E ]Pao/'ae
+ a p../-.. + 2 0(1. - p../n..)'°‘£. <3-37a>
1.. = E p..I-.. - I E + 2.0 1.1,)... - .,
(3.37b)

where E = 3TANH(P1./I1.’p;1/flxg)l(ng/lgg’pgg/Igg)

Comparing the terms in the brackets with the equivalent
parameters of these R elements defined by the chordal approximation it
turns out that they are identical. At a certain global time the state
equations in matrix form would be

i)... -(a,+a,+a.) a, Us... 0 p... r.-p,./a,1‘-‘) 0 f, (3.33.)
it; . R. -(R,+R.) 0 l/m1 p11 0 l e, (3.38b)
[A]

where R1. R,. R,. R4 are defined by Eqs.(3.34a.b.c.d).

With the equivalent linear resistances. an estimate of the
instantaneous dynamics of the nonlinear system can be obtained by the

extraction of the eigenvalues from the above [A] matrix. The [R]

67

matrix is analogous to the one of linear case described earlier.
Although this analogous [R] matrix can not be used for direct simula-
tion. it may serve to predict the local dynamics in the two-time-scale
integration technique. It will help in parameter selection for the
dynamicizers and the local integration controls. Therefore. such an

[R] matrix is named the reference [R] matrix.

By using the instantaneous values for the bond variables or state
variables available from the previous global time step. a set of
instantaneous equivalent linear resistances can be computed. In this
example. at a certain global time step the state variables p1. and pu
‘10 known 38 Pa. and 511. The reference [R] matrix would be calculat-

ed as follows

[R] = - [(T,-511/511)°°’ + in];u + E ] E
E -[E+2.0]
where T.. 510.511. 3,. and 3,, are the values at the previous global

time step. and E is as in Eq.(3.37).

In many cases a set of R parameters could be used for several
global time steps before re-evaluation would become necessary. This
must be judged from the rate at which the local field input variables

are changing.

At steady state the equivalent resistance parameters virtually

may be the functions of only the input vector to the R-field; namely

68

R. = 1.0)) (3.39)

If the change of the input vector is not too large the equivalent
resistance Ri will not change dramatically. provided the nonlinear
dissipative element is well-behaved; for instance. the effort is a low
order power of flow. For many physical systems. the nonlinear dissi-
pative elements possess fairly moderate characteristics. The changes
in equivalent parameters ever the dynamic range of interest are not
very large. In the mechanical systems such nonlinear dissipation
effects come from static friction. columb friction and other nonlinear
frictions which have small degree of nonlinearity. Therefore. the
reference [R] matrix will not change the local dynamics too much and

will not require the freqent re-selection of added dynamicizers.

The process described in the preceding paragraphs can be repeated
as the global time variable increases; that is. the free parameters
could be reselected intermittently throughout the global simulation in

prescribed fashion.

The technique of parameter selection for the n-th-order nonlinear

subsystem is summarized below:

1) At t=t.. set all free parameters to be an arbitrary constant.
say. unity. set arbitrary initial conditions for the added dynamic
elements and choose a small time step and a long duration for local

time scale. then integrate the subsystem to steady state.

69

2) Compute the instantaneous equivalent linear resistances for

each nonlinear dissipative element in the augmented subsystem.

3) Form the reference [R] or [0] matrix.

4) Select the optimal free parameters following the general rule.

5) Compute the eigenvalues of -[G][R] or -[R][N“] matrix. reset

the local time scale.

7) Integrate the subsystem to the steady state in the new local

time scale.

8) Compute the output vector of the subsystem.

9) Continue with the global integration forward one step.

10) Repeat the parameter selection process according to the error

control criteria.

70

4 NUMERICAL EXAMPLE

To illustrate the effectiveness of the bond graph augmentation
method discussed in the preceding chapters. a demonstration example is
presented in Figure 4-1. The subroutine PART identified the system
and partitioned it into one dynamic field and one static field. Then
subroutine AUGMEN automatically determined which nodes should be aug-
mented and what kind of dynamicizers should be added. After the
augmentation process was finished. the augmentation information was
printed out for inspection (Figure 4-2). The augmented subsystem and
the dynamic subsystem are shown in Figure 4-3. The bonds and the
nodes in both subsystems have been renumbered regarding the initial

descending numbering sequence.

For different purposes a linear case and a nonlinear case have
been simulated. The existing program for linear systems can be used

as an examiner.
4.1 Linear Case

If the all dissipative elements in the erield are linear. then

the state equations of the augmented subsystem are

i)... -(a,+a.) a, 1n... 0 p... 0 -1 r. (4.1)
..-

1311 R, -(R1+R,+R,) O 1/mu p11 R1 0 e,

71

(8) (9) (10) (11) (12) (13) (14) (7)

SPF—1‘0 —-9 1———1° 0—-11 1|-—-12 0—-|13 1——-|7 I

R R R R C R
(1) (2) (3) (4) (5) (6)

Figure 4-1 The global bond graph

model for the example

72

IL HE ARE 1 R~FIELDS IN THE GIVEN SYSTEM

II FIELD NUI'IUI.R 1
NIIIJILS
9 O B 1 9
8 SF 8
I R 1
IO 1 9 2 10
a R 2
11 O 10 .3 11
3 R 3
12 1 11 4 12
4 R 4
13 SE 12
III)“ [S
12

THERE ARE 0 JUNCTION STRUCTURE COMPLEXES
THERE ARE 1 DYNAMIC FIELDS IN THE GIVEN SYSTEM

DYNAMIC FIELD NUMBER 1
I'J‘JDES ,
13 O 12 5 13
12 SF 12
5 C 5
l4 1 13 6 7
6 R 6
7 I 7
PORTS
12

ADD C/I ELEMENT TO THE ACTING JUNCTIONS

ADD 1 ELEMENT TO NEH NODE B (THE OLD NODE 12 )
ADD C/I ELEMENT TO THE SINGLE-TYPE JUNCTINS

ADD I ELEMENT TO NEW MODE 4 (THE OLD NODE 10 )

THERE ARE 2 STATE VARIABLES
AND .2 INPUT VARIABLES.

THE STATE VECTOR...
x 1 = P 10
x 2 =~P 11

THE INPUT VECTOR...
u 1 = F 5
u 2 = E 9

Figure 4-2 The Augmentation Result for the Example

73

P f
10 5
X‘ 3 Us 8 v: 2 f9
1’11 ’9
( s )

[p3]

SFI'L-OL‘Tl-i‘ll :m
I ’T
. C :k R

[91]
(b)

Figure 4-3 The subsystems
s) Augmented subsystem
b) Dynamic subsystem

74

If all the elements in the dynamic field are linear. the state

equations of the dynamic sybsystem are

as = fa '- 93,”: (4.23)

fit ‘ qul ’ Rspa/”a (4.2b)

Assume the parameters of the elements are

R1 = R: = 10.0 a, - n: - 20.0
R, - n: = 5.0 n. - n: - 10.0
k, = 12 = 20.0 n. = 32 = 5.0

m, I m? = 2.0
where the superscripts s and d denote the augmented subsystem and the

dynamic subsystem. respectively. The [3] matrix of the augmented subsystem is

[R] = -15.0 5.0

5.0 -35.0

According to the general rule ve choose 15 and 35 for parameters

l1. and ‘11: respectively. so the [A] matrix appears

[A] ' -1.0 0.1429
0.3333 -1.0
The eigenvalues of the [A] matrix are
a, - -O.7818
a, = -1.2182

We may set the local time scale as follows

75

At = 1"3‘a1n’ - l/(3‘1.2182) - 0.273 sec

tf - 5/s_.x - 5/0.7818 - 6.395 sec

If you want to change the scale by a factor n. you may change the
parameters l3. and mg; by the factor. For instance. let m1. and m1,

be 1.5 and 3.5. respectively. then

‘1 = -70818
‘3 x -12.182
then the time interval and integral duration will be

At = 0.0275 sec

tf = 0.6395 sec

In the execution we set At = 0.025 sec and tf = 0.65 sec.

The global time scale may also be determined from the eigenstruc-

ture. 81,, = -1.95 t 3.114j. as follows

AT = 0.05 - 0.10 sec

Tf = 5.0 sec

Suppose that all the initial conditions are equal to zero at T =
O. and the flow input from source SF8 is a constant. say. f: , 10,0,

Since T = 0. q: = 0. and p, - 0. therefore e: = e? . 0 . These
two subsystems interact through bond 12 in the initial system

model.DOhat is the output of the augmented subsystem. f:, is the input

76

t0 th° dynamic subsystem. ‘2. while the output of the latter. of. is
the input to the former, 0:. Besides the input to a augmented subsys-
tem at present global time step from the adjacent dynamic subsystem
would be the output of the latter at the just previous step. 80. at
global time T 8 0.05 sec. the input e: s 0. Integrating in the local

time scale that has been set earlier yields the output f: at T = 0.05.

f: 3 Sic/nio° (4.3)

It will be the input to the dynamic subsystem at T = 0.05.
Figure 4-4 shows the local time scale integration process at the glo-
bal time T - 0.05. and Figure 4-5 displays the next local integration.
The dynamic subsystem then integrates one time step in the global time
scale. The steady states of the augmented subsystem serve as initial
conditions for next local-time-scale integration and also the states
of dynamic subsystems become the initial condition for the next
global-time-scale integration. These two subsystems interact in this
manner throughout the entire global time so that the system simulation

can be realized.

To demonstrate that this procedure can lead to correct results.
we have used the DIFFEQI12] package that can solve either linear and
nonlinear differential equations. The numerical results and a plot
are shown in Figure 4-6 and Figure 4-7. respectively. Comparing these
with the results obtained by the ENPORT‘SI13] package (Figure 4-8 and

Figure 4-9) shows that they are almost identical. The small differ-

77

con 06‘»

m

uuuuu emu uo>o momuuuuoumm alumna-cod vuv nonumm

N

noun alumna-moan

 

 

«an

. 3..

 

 

m.N

o.m

mK

9.0—

m.N_

78

oo- oluu m

emcee. on» nose coma-uuoum« alumna-ecu n1. cunnmu

N

noun oluutuaacuu

_

 

—

-

 

2..

d.“

 

9.0

5.9.

m.N_

”Huuuuuuuuwuuuoommqqo

99999999999999000999999999900000o
()flmVO‘UIbUI'Jl'Jflt-‘O .

TIME

.SODCE-Ol
.lCOCE+OO
.ISOCE+OO
.2000E+OO

SOOE+OO

Iaoocs+oo

SOOE+OO

OOE+01

0cm)
rmnm
+++
0cm:
HHH

on
mnl
++
00
HH

OE+Ol
CE+01

00
mm
+
O
..a

08000000000
m '

+ +

o o

eds-e

OE+01

QCKMDOCKMDOCKMDOCKMDOCKMDOCKNDOCKMDOCKMDO

P10

.49BOE-01
.9497E-01
.1349E+OO
.1693E+OO
.19806400
.2212E+OO
.2391E+OO
.25218+OO
.ZbObE+OO
.26528+OO
.26658+OO
.2650E+OO
.26128+OO
.25546+OO
.24B9E+OO
.24138+OO
.233SE+OO
.22525+OO
.217BE+OO
.2098E+OO
.2029E+OO
.1966E+OO
.l9llE+OO
.leAE+OO
.1825E4OO
.1794E+OO
.1752E+OO
.l736E+OO
.1739E400
.1754E+OO
.177bE+OO
.lBOOE+OO

1822E+OO

79

OCKMDOCKMDOCKNDOCKMDOCKMDOCKMDOCKMDOCKMDO

P11

.2394E-01
.89745-01
.187BE+OO
.3093E+OO
.4462E+OO
.59llE+OO
.7384E+OO
.8829E+OO
.1020E+01
.ll485r01
.1263g+01
.1365E+Ol
.1452E+Ol
.1524E+01
.15825+Ol
.1627E+01
.16598+Ol
.1680E+Ol
.lb9OE+Ol
.169ZE+01
.1687E+Ol
.16766+Ol
.1661E+Ol
.16438+Ol
.1623E+01
.1602E+Ol
.lelE+Ol
.1524E+Ol
.1494£+Ol
.1473E+01
.145?E+01
.1453E+Ol
.1452E+01

IOCKMDOCEMDOCKflDOCKMDOCHMDOCKMDOCKMDOCK”?O

Figure 4-6 Numerical result of linear case

E9

.99605+OD
.1899E+Ol

269BE+01

.3386£+Ol
.39blE+Ol
.4424E+Ol
.47825+01
.50428+Ol
.52138+Ol
.5305E+01
,5329E+Ol
.52?9E+Ol
.5224E+Ol
.5108E+01
.4979E+Ol
.4827E+Ol
.4566E+Ol

504E+Ol
4346E+Ol

.4196E+Ol
.4057E+Ol
.39335+Ol
.8822E*Ol
.3728E+Ol
.36505+Ol
.3588E+Ol

504E+Ol

.3471E+Ol
.3478E+Ol

35085+Ol

Iass1s+01

3599E+Ol

.3644E+01

obtained

by the two-time-scale integration

80

mcU-unoumm cane-leluulcsu 23 ha coma-«no oomonuou 3|.th TV 93.:

 

 

 

 

 

 

 

 

 

 

99999999999999999999999999999999999999990

TIME

.OOOE

00

01

OCKMDOCKMDOCKMDOCI
HHMHHHHHHHHHHH

999999999999999999999popppooopooooooooooo

P10

.OOOOOE
.480955-01
.918625-01
.13073E
.16437E
.19269E
.21577E
.23382E
.24718E
.25629E
.2616CE

26364E

.26292E
.259945
.25521E

81

00

P11
OOOOOE

00

233865-01

13411E

.14294E
.15036E
.1564OE
.16112E
.16461E
.16697E

O.
0.
0
0
O
0
O
0
O
0
O
0
O
0
0
O
0
O
0
O.
O.
O.
0
O
0
0
O
O
O
O
O
0
O
O
O
0
0
O
0
0
0

497E

VQQOHUUVQHU

.87213E-01
.18237E
.300385
.43361E
.57527E
.719545
.861555
.99736E
.11239E
.123905

Figure 4-8 Numerical result obtained by the ENFORTBS

82

ss.m

«Taxonzu emu hm conquumc o-moauou cum-ch: «Iv shaman

......zmx .llldx

 

 

 

s

\1

Q
Q

”ozmmmn

 

s¢.s

 

sm.s

 

sm.~

 

 

 

 

 

 

 

 

 

5 mad» w4<um wth
am." am." sm.s sv.s
......-........4..........
...4.......1.......

 

 

so.“
5 mad“
mm<um »

ss.m

 

 

83

ence (about 1.0 - 2.0 i ) comes from the different algorithms employed
in these two programs. It. however. has clearly illustrated how the

two-time-scale integration technique works and how well.

4.2 Nonlinear Case

If some of the dissipative elements in the augmented subsystem
are nonlinear. the simulation process will follow the steps summarized
in the preceding chapter. Suppose the dissipation functions are the

same as Eqs.(3.34a.b.c.d). The state equations are rewritten as below

01. = 3TANH(p11/m.,-p,./m1.) - 2.0p,./m1° - e, (4.5a)
P11 ‘ 2'OSIGN(f5—p11/n11)(fl-p11/m11)1.a ' SIGN(Pie/lie)(Pie/‘ie)3
- 3TANB(p11/I11-p1./ll1.)

(4.5b)

The reference [R] matrix may be

-(R,+R.) 3:
[R] . -[ E + 2.0 1 a
E -l<f.-p../-..)°-'+p../-..+ E 1
R. -(n,+R,+R,)

where E g BT‘NH(p11/'11 ' Pie/“io)/(Paa/“xi ’ Pie/“10)

The initial partmetor' I... m,, and the initial conditions p...

pu may be any reasonable numbers as long as there is no zero divide

implied. say. m1° = 1.0. 11111 8 1.0. p,. 8 1.0 and p11 8 2.0.

84

Choose At = 0.01 and tf I 3.0 as the local time scale.
Integrating in the local time scale produced the steady state values
5,, and S... and then all the bond variables can be calculated from
the constraint equations. They are used to compute the instantaneous
equivalent linear resistances for R1. 3,. R, and 3.. The referece [R]

matrix. then is formed as below

[R] I -23.351 0.7223
0.7223 -2.7223

For a smaller time scale. we set

I;. I 2.3351
m“ I 0.27223
Then the [A] matrix is
[A] I -10.0 2.653
0.3093 -10.0
The eigenvalues of [A] are
s1 I -9.0941
'3 I —10.9659
Reset the local time scale by
At I 1/(3‘10.9059) I 0.03056 I) 0.030 sec
tf I 5/9.0941 I 0.5498 I) 0.60 sec
For comparison the eigenvalues and the local-time-scale have been
computed when parameters m1. and m,; are arbitrarily set
to be one. they are listed below:
I, I -2.697
sa I -23.376

and At I 0.01425 sec

85

tf I 1.854 sec

Figure 4-10 displays the local integration process when the free
parameters are selected arbitrarily. while Figure 4-11 shows the local
integration at the same global time after the parameters have been
optimized. From these plots we can see that the optimal parameters
make the "diffusion” process quicker and the local integration steps
fewer. The remaining process would be the same as that in the linear
case discribed in preceding section. The only difference is that the
nonlinear case needs to control the local time scale by detecting the
change rate of the subsystem input vector or the equivalent resis-
tances. Fortunately. the equivalent resistances. R1, 2,, R, and 3‘ in
the particular example change only a little. so the local time scale
needs not to be modified in the later 20 global time steps. The
numerical output and a plot of this nonlinear case are shown in Figure

4-12 and Figure 4-13. respectively.

Comparing Figure 4-13 with Figure 4-9 obtained in linear case
shows that they have a similar appearance aside from the different
y-axis scales. It means that both cases have almost the same dynamic
characteristics. This may be explained by the fact that the nonlinear
dissipative field changes are very small in the particular system

because of the constant input SF8.

86

annual-nan huuuuqmua duds momunuuoumm "coca onlv anemia

no. alga

 

2a

 

3..

 

87

 

auouolauaa «engage mums memuuumoumm mecca «21¢ shaman

con old»

 

 

3a

 

Sq

 

ococcoooocooooocccoo

Fja-He-Ie—He-Ia-se-e

IIME

.lOOOE+00

EOOOE+OO
SOGOE+OO
4000E+OO
SOOOE+OO
6000E+OO

.7000E+OO

BOOOE+OO
9000E+OO
IOCOE+01
llCOE+Ol

rrrn
++
00
e-IH

+01

o~cmv0~mmam
000000000
000000000

++++++
oooooo
Hus—oped»-

Figure 4-12 Numerical result obtained by the

oooogoogpppooooooooo

P10

.8745E-01
.8993E‘01
.BS75E-01
.79OBE~OI
.7582E-01
.72115-01
.7364E-01
.7110E-Ol
.704OE-01

6993E-01
695BE~01
6936E-01
6924E-01

.6909E-0l
.6906E-01

6903E-01

.6897E-01
.6899E-01
.6895E-01
.6895E'01

two-time-soale integration

.8
.22135+OO
.3266E+OO
.398BE+00
.4479E+OO
.50385+OO

00000000000000000000

P11
121E"Ol

4812E+OO

.5191E+OO

5295E+OO

.53665+OO
.5414E+oo
.544bE+00
.5468E+00
.5482E+OO
.54925+00
.5499E+00
.5504E+00
.5507E+oo
.55095+oo
.5509E+oo

99999999999999909999

HHHHMHMHHs—wubsb-kus—suw
UUUUUUUUL)UU-bb&bmmo~q
uvmummmmmeoomunumuo
oaoowummumommumomm
mmmmmmmummmmmmmmm
+++++++++++++++++
OOOOOOOOOOOOOOOOO
WHHHHHHHHHHHHHHHH

89

once woomuunon am» no canon-cu oglumha mnlv onen«m

 

 

 

 

 

 

 

90

5 CONCLUSION

5.1 Summary

The bond graph augmentation method for simulation of a dynamic
system containing coupled nonlinear dissipative fields has been

developed.

The causality assignment process is very useful for identifying
Refields. The order of augmented subsystems has been briefly dis-
cussed and it has been shown that to get the least augmentation order
and to make the parameter selection as easy as possible. adding only
one type of dynamic element ( C or I ) is desirable. An augmentation
sequence has been suggested. By adding only one type of dynamic ele-
ment the augmented subsystem will reach its steady state in a
"diffusion” manner. Tb obtain the best computation efficiency the
selection of parameters for the dynamicizers is the key. The discus-
sion on parameter selection in linear subsystems has been extended to
nonlinear subsystems through the concepts of instantaneous equivalent
linear resistance and the reference [R] matrix. A general rule for
parameter selection for arbitrary n-th-order subsystems has been sug-
gested and numerically verified by several cases with different orders

and numerical conditions.

The partition of dynamic fields and dissipative fields and the

augmentation process are automatically accomplished by the subroutines

91

PART and AUGMBT.

An example with linear and nonlinear cases has demonstrated the
effectiveness of the new approach. The results obtained by this
approach and by the existing ENFORT-S program from linear case have
been compared. It shows that the two-time-scale integration technique
is valid for dynamic study. As the dissipative effects in the same
system are reset to be nonlinear. the simulation result has shown the

expected system dynamic behavior.

The approach to simulate a system with coupled nonlinear dissipa-

tive fields proposed above has several advantages. Among them are:

l) The parallelism in this approach is well-suited for parallel

computers.

2) The manner in which we break up the whole problem into parts

may be used to simulate large systems in a small capacity computer.

3) The coupled algebraic loop in system equations will be solved
and the system dynamics can be investigated without necessarily

approximating the nonlinear dissipative effect as a linear one.

On the other hand. the use of this approach based on only C or

only I augmentation caries with it certain possible limitations:

92

1) For the adverse numerical condition of a dissipative field(its
"spread” of eigenvalues is close to one). the computational efficiency

would deteriorate.

2) The general rule for parameter selection may not apply to cer-
tain nonlinear static fields containing junction lo0ps or GY elements
because sometimes it is impossible to complete the causality by adding

only C or I elements.

5.2 Future DevelOpment

The complete implementation of the approach is not possible until
the program NONLIN. which can simulate uncoupled nonlinear systems. is
available for use. This program is currently under development. The
remaining work mainly is to design the control of the data flow among

subsystems and the data assembly.

Basically the mixed C and I augmentation is not preferred becuse
it might cause complex eigenvalues. But it still remains for futher

investigation.

Some interesting problems are the sensitivity of the eigenstruc-
ture of an augmented subsystem to the change of its input vector and
the influence of the augmentation order .on the sensitivity.
Investigation into these problems may benefit the control of the local

time scales and therefore computational efficiency.

93

Until this stage. it is still too early to give a definite con-
clusion on the computational efficiency of this approach. The
efficiency is commonly defined to be inversely proportional to the CPU
time consummed during the execution of the computer program. One
needs to compare the CPU time that is used to obtain the solution for
different class problems by this approach with those used by other
traditional methods. such as mentioned in chapter 1. It would be pos-

sible when the program NONLIN is available.

REFERENCES

94

1). CSMP III . lhite Plains. New Yoyk. IBM. 1975.

2). R. N. Nilsen. CSSL-IV . Chatsworth. CA.. Simulation Service.
1977.

3). G. A. Korn. J. V. Iait. Digital Continuous System Simulation .
Bnglewood Cliff. New Jersey. Prentice-Hall. 1978.

4). J. C. Bowers and S. R. Sedore. SCEPTRE . Englewood Cliff.
Prentice-Hall. 1971.

5). G. D. Byrne and C. A. Ball. Numerical Solutions of Systems of
Nonlinear Algebraic Equations . Academic Press. New York. 1973.

6). C. Wolfe . Phillip. The Secant lethod for Simultaneous Nonlinear
Equations . Communications of the AC“. v01. 2. 1959.

7). K. I. Brown. A Quadratically Convergent Newton-like lethod Based
Upon Gaussian Elimination . SIAM Journal on Numerical Analysis.
Vol.6. 4. 1969.

8). The INSL Library . International lathematical & Statistical
Libraries. Houston. Texas. 1982.

9). R. C. Rosenberg. D. C. KarnOpp. Introduction to Physical System
Dynamics . New York. NcGraw-Uill. 1983.

10). P. L. Graf. Dynamic Augmentation of Dissipation Algebraic
Loops . I.S. Thesis. Michigan State University. 1981.

11). A. R. Gourlay and G. A. Iatson. Computational lethods for
latrix Eigenproblems . New York. John Wiley & Sons. 1973.

12). Case Center User's Guide . A. H. Case Center for CAD/CAM.

95

Richigan State University. 1982.
13). R. C. Rosenberg. ENFORT-S User's lanual . Iichigan State

University. 1983.

APPENDIX A

96

(ENFORTG)

(CAUSE)

FAR AUGMEN

| +|
l l l I I] l

RFIELD DFIELD PRDFID REARRG UNINAI BAKNAM ADELMT

 

 

CALLING TREE

APPENDIX B

97

COPART‘...0......Ot......OOOOOOOOOOOOOOO............OOOOOOOOOOOOOO

C
SUBROUTINE PART

C

C -- PROGRAMMER -- TONG ZUOU. AUG. 1983

C

C -- PART FINDS AND SEPARATES TEE RIFIELDS AND THE JUNCTION
C STRUCTURE COMPLEXES FROM TRE.ORIGINAL BOND-GRAPH.
C

C MR -- RPFIELD INDEX

C MD -- DYNAMIC FIELD INDEX.

C NJ --- CAUSAL BRANCH INDEX

C JX -- JUCTION STRUCTURE COMPLEX INDEX

C

can DECLARATIONS
c
smssx'r syncsx
$INSERT CAUSBK
smsam UTILBK
$INSERT mm:
c
C03000OOOOOOOOOOOOOOOOOO00.00.0000.000000COOOOOOOOOOOO000000.00...
c

113:0

1m=0

NJIO

1301=0

IBR1=0

NELR=O

NELD=0

INN=0

NBDZINBD‘Z

DO 2 II1.NBD2
ICMXT(I)=ICMX(I)
CONTINUE

2
C
C --- FIND THE FIRST JUNCTION
C
5

DO 10 N1I1.NEL
IF(IELLST(N1).GE.6) GOTO 15
10 CONTINUE

15 NPlINPTR(N1)
NP2INPTR(N1+1)-1
C
C --- DETERMINE IF ALL THE BONDS INCIDENT TO THE JUNCTION
C ARE CAUSAL ASSIGNED
C
DO 25 IINP1.NP2

IF(ICMXT(I).EO.1) GOTO 30
25 CONTINUE

c-..
C

100

110

C

98

GET THE ALL INFORMATION FOR THE PRESENT DYNAMIC FIELD

IF (MD.EQ.0) MD=1
CALL DFIELD(MD.MR.N1.NP1.NP2)
GOTO 35

GET THE ALL INFORMATION FOR THE PRESENT RPFIELD

MR=MR+1
CALL RFIELD (MD.MR.N1.NP1.NP2)

EXTEND TO THE NEXT ADJUCENT FIELD

IF(IXBD.EQ.0) GOTO 40
N1IJBD(IXBD)
IXBDIIXBD-l

GOTO 15

IF (IXBR.E0.0) GOTO 90
N1IJBR(IXBR)
IXBRIIXBR-I

GOTO 15

FIND THE JUNCTION STRUCTURE COMPLEX FROM R-FIELDS

JX=0

DO 110 I=1.MR

LILR(I)

DO 100 JI1.L

IF (IRLST(I.J).EQ.3) GOTO 110
CONTINUE

JX=JX+I

JCOMX(IX)=I

CONTINUE

INOIINN

C--REARANGE THE DATA ARRAYS

C

120

D0 120 II1.MR

MDlIO

IBRIIIBR(I)

NELRILR(I)

CALL REARRG ( I.MD1.IBR1.NELR.IBMXR.NBIMXR.IBMXRN.

+ NBIXRN.INDXRE)

CONTINUE

DO 130 II1.MD

MRIIO

IBDIIIBD(I)

NELDILD(I)

CALL REARRG ( MR1.I.IBD1.NELD.IBMXD.NBIMXD.IBMXDN.

+ NBIXDN.INDXDE)

99

130 CONTINUE

C
C -- PRINT R-FIELDS AND JUNCTION COMPLEXES
C
CALL PRRFLD(MR.MD)

C

MRTIMR

MDTIMD

RETURN

END
C
C
C
CRFIELDsssessasseessseseesessseesesssesseesseseeesssessessssseee
C

SUBROUTINE RFIELD (MD.MR.N1.NP1.NP2)

-- RFIELD GETS ALL INFORMATION FOR THE RPFIELDS( zUNCTION
STRUCTURE COMPLEX IS TREATED AS A R-FIELD ) FROM THE ORIGINAL
BOND-GRAPH.

INPUT --- N1.NP1.NP2

OF THE BEGEINING JUNCTION OF THE RrFIELD
OUTPUT -- N1.NP1.NP2

OF THE END JUNCTION OF THE R-FIELD

L -- NODE INDEX

IB -- BOND INDEX

IP -- PORT INDEX

NUMJ -- BRANCH INDEX

INDXRE 4- NODE SERIEL NUMBER

IRLST -- ELEMENT TYPE LIST

IRNAM -- ELEMENT NAM LIST

NPTRR -- POINTER LIST FOR NBIMXR(START OF BOND GROUP)
NBIMXR --- LIST OF BONDS INCIDENT'TO EACH ELEMENT
JUNCTR --- LIST OF BRANCHES IN THE R-FIELD

LR -- NODE NUMBER STORRAGE ARRAY FOR EACH RrFIELD
IBR -- BOND NUMBER STORAGE ARRAY

IPR -- PORT NUMBER STORAGE ARRAY

fiffifbf)(ififlfhfififlfifif5f5fi(Bflfifhfifififi

cu. DECLARATIONS
c
$INSERT smear
smsan'r CAUSBX
$INSERT UTILBK
$INSERT PARTBK
c -
c
c0000...OOOOOOOOOOOOO0.000000000000000....OOOOOOOOCOOOOOOOOOOOOO
c
c -- INITIALIZATION
c
L-o

55
C

100

IBIO

IPIO
NUMJIO
NN2=0
NBD2=NBD‘2

GET JUNCTION DATA

L=L+1

NPTRR<MR.L)IIB+1

IIIIB

DO 20 I=NP1.NP2

IF (ICMXT(I).E0.0) RETURN
IBIII+I-NP1+1
NBIMXR(MR.IB)INBIMX(I)
ICMXR(MR.IB)=ICMX(I)
IBOND=NBIMX(I)
IBMXR(MR.IBOND.1)IIBMX(IBOND.1)
IBMXR(MR.IBOND.2)IIBMX(IBOND.2)
CONTINUE.

IRLST(MR.L)‘IELLST(N1)
IRLNAM(MR.L)=IELNAM(N1)
IF(L.EQ.1) GOTO 50

DO 30 I=1.LP1
IF(INDXRE(MR.I).EO.INDXRE(MR.L)) GOTO 40
CONTINUE

GOTO 50

IBIIB-NP2+NP1

LIL-1

GOTO 185

GET ADJOINING NODE DATA

DO 180 IINP1.NP2
IBONDINBIMX(I)
NAIIBMX(IBOND.1)
NBIIBMX(IBOND.2)
N2INA

IF(N1.EQ.NA) N2=NB

NUM=0

D0 55 JIl.60
IF(NBIMXR(MR.J).NE.IBOND) GOTO 55
NUM=NUM+1

IF(NUM.NE.2) GOTO 55

ICMXT(I)I0

GOTO 180

CONTINUE

C --- SAVE BOUNDARY INFORMATION

C

101

IF (IELLST(N2).EQ.3) GOTO 95
IF (IELLST(N2).EQ.4) GOTO 95
IF (IELLST(N2).EQ.5) GOTO 95
IF (ICNXT(I).EQ.1) GOTO 95

NN2=NN2+1

ID=ND+1

DO 60 J=1.NBD2

IF(J.EQ.I) GOTO 60
IF(NBIMX(J).E0.IBOND) GOTO 65
CONTINUE
IF(IELLST(N2).EQ.1.0R.IELLST(N2).E0.2) ICNXT(J)=O
IF(ICNXT(I).EQ.6) GOTO 70
IF(ICXXT(I).EQ.3) GOTO 75

SINGLE C- ELEMENT

ICMXT(I)-O

ISOURL(NN2.1)-5
ISOURL(NN2.2)=4
ISOURN(NN2.1)='SF'
ISOURN(NN2.2)-'SE'
Icuxn(Mn,1)-6

Icuxn(un.2)=3
IF(IELLST(N2).GT.5) GOTO 85
GOTO 80

C-- SINGLE I- ELEMENT

C
75

ICMXT(I)=O

ISOURL(NN2.1)'4
ISOURL(NN2.2)=5
ISOURN(NN2.1)='SE'
ISOURN(NN2.2)-'SF'
ICNXD(MD.1)=3

ICIXD(ND.2)=6
IF(IELLST(N2).GT.5) GOTO 85

GET THE DATA FOR THE SINGLE DYNAMIC SUBSYSTEH

NPTRD(MD.1)=1
NPTRD(MD.2)=2
NPTRD(MD.3)=3
INDXDE(ND.1)=N1
INDXDE(MD.2)=N2
IDLST(ID.2)-IELLST(N2)
IDNAM(ND.2)=IELNAM(N2)
IDLST(MD.1)RISOURL(NN2.2)
IDNAM(MD.1)=ISOURN(NN2.2)
NBIMXD(ND.1)'NBIMX(I)
NBIMXD(ND.2)=NBIHX(I)
IBNXD(ID.1.1)-N1
IBNXD(ID,1.2)'N2

c ......

85

90

100
105

110

115

102

IBD(ID)=2
IBD1=2
LD(ND)=2
NELD=2

CHECK IF THE JUNCTION HAS BEEN TAKEN

IF(IELLST(N2).EQ.1.0R.IELLST(N2).EQ.2) GOTO 90
IF(ICNXT(J).EQ.0) ID8lD-1
NUNJ=NUNJ+1
JUNCTR(NUNJ)'N2

IP-IP+1
IPORTR(NR.IP)*IBOND

L=L+1

IB=IB+1

NPI'RR(MR,L)=IB
INDXRE(NR.L)=N2
IRLST(NR.L)=ISOURL(NN2.1)
IRLNAM(MR.L)=ISOURN(NN2.1)
NBIMXR(NR.IB)=NBINX(I)
ICNXR(MR.IB)=ICNX(J)
INN=INN+1
INTER(INN.1)=IBOND
INTER(INN.2)=IBOND
INTER(INN.3)'IR
INTER(INN.4)=ND
INTER(INN.5)=ISOURL(NN2.1)
GOTO 180

ICMXT(I)=0
IF(IELLST(N2).GT.5) GOTO 110
IF(IELLST(N2).EQ.1.0R.IELLST(N2).EO.2) GOTO 180

L=L+1

IB=IB+1

NPTRR(MR.L)=IB
INDIRE(HR.L)=N2
IRLST(MR.L)-IELLST(N2)
IELNAM(NR.L)=IELNAN(N2)
NBIMXR(IR.IB)=NBIMX(I)

DO 100 JSI.NBD2

IF (J.EQ.I) GOTO 100

IF (NBIMX(J).EQ.IBOND) GOTO 105
CONTINUE
ICIXR(NR.IB)=ICIX(J)
ICMXT(J)=0
IF(IELLST(N2).LT.6) GOTO 180

D0 115 K=1.L

IF (INDXRE(NR.K).EQ.N2) GOTO 180
CONTINUE
NUNJ=NUNJ+1

103

JUNCTR (NUNJ)-=N2
180 CONTINUE
NT=N1
C
C --— NEXT JUNCTION
C

200 IF(NUHJ.EQ.0) GOTO 210
185 N1=JUNCTR(NUMJ)
NUNJ=NUMJ-1
NP1=NPTR(N1)
NP2=NPTR(N1+1)-1
DO 205 I=NP1.NP2
IF(ICMXT(I).EQ.1) GOTO 10
IF(ICMXT(I).EQ.0) GOTO 200
205 CONTINUE

IXBD=IXBD+1
JBD(IXBD)=N1
GOTO 200
C
C -- STORE THE SIZE DATA
C
210 IPR(NR)=IP
LR(NR)=L
IBR(NR)=IB
IBR1=IB
NELR=L
NPTRR(MR.L+1)=IB+1
C
C
RETURN
END
C
C
C
copRRanoot0.0.0.0000c.0000oooooo¢¢ocoooccacoooooo¢¢ocotto0000000
C

SUBROUTINE PRRFLD(NR.ND)
C
C -- PRRFLD PRINTS OUT THE INFIRNATION ABOUT THE RFFIELDS AND
C JUNCTION STRUCTURE COMPLEXES.

C

C -- DECLARATION
C

SINSERT PARTBX

C

COOOOOOOOOOOOOOOOOOOOOOOO.........OOOOOCOOOOO0.0.0.0...00.0.0...

C
C -- PRINT RFFIELD
C
OPEN(UNIT=5.FILE='SBOW')
C
IR=NR-JX
'RITE(‘.1000) IR

1000

950

900

960

1010

1020

150
5030

1050

1060
100

c_..._

300

1100

1210

104

FORKAT(/'THERE ARE '.IZ.' R-FIELDS IN THE GIVEN SYSTEH')
IF(IR.EQ.O) GOTO 300

DO 100 I‘1.NR

IF<JX.NE.O) GOTO 950

N80

GOTO 960

D0 900 N‘I.JX

IF(I.NE.JCONX(N)) GOTO 900

GOTO 100

CONTINUE

IN=I'N

'RITE(‘:IOIO) IN

FORMAT(/'ErFIELD NUNBER '.I4.' :')
'RITE(‘.1020)

FORMAT(/'NOOES :')

DO 150 K‘1.LR(IN)

II'NPTRRU .K)

I2=NPTRR(I.K+1)-1

'RITE(‘.5030) INDXRE(I.K).IRLNAN(I.K).(NBINXR(I.J).J=I1.I2)
FORMAT(7X.I3.21.A4.ZX.5I4)
'RITE(..IOSO)

FORMAT(/.'PORTS :'ol)

'31TE(‘.1060) (IPORTR(I.J).J=I.IPR(I))
FORMAT<7xaI3)

CONTINUE

 

PRINT JUNCTION STRUCTURE COMPLEXES

CONTINUE

'RITE(‘.IIOO)JX

FORMAT(/'THERE ARE '.I2.' JUNCTION STRUCTURE CONPLEXES')
IF(JX.EQ.O) GOTO 260

DO 200 I‘lpJX

JXN=JCOMX(I)

'EITE(..1210) I

FORMAT(/'THE JUNCTION STRUCTURE COUPLE! NUMBER 'aIZ.’ :')
'RITE(..IOZO)

DO 250 K‘I.LE(JXN)

II'NPTRR(I.K)

IZ'NPTRR(I.K+I)’1

'RITE(‘.5030) INDIRE(I.K).IRLNAN(I.K).(NBINXR(I.J).J=I1.IZ)
'RITE(‘.IOSO)

'RITE(‘.IOOO)(IPORTEIJXN.J).J‘I.IPR(I))

CONTINUE

PRINT DYNAMIC FIELDS

IRITE(‘.1300) ID

FORMAT(/'TEERE ARE '.12.' DYNAMIC FIELDS IN THE GIVEN SYSTEM')
DO 400 I=1.ND

'RITE(‘.1310) I

FORMAT(/'DYNAMIC FIELD NUMBER '.I4.' :')

'RITE(‘.1020)

105

D0 310 R=1.LD(I)
II-NP'I'RD(I.X)
IZ-NPTRD(I.K+1)-1
310 'RITE(’.5030) INDXDE(I.K).IDNAM(I.X).(NBIMXD(I.J).J=II.IZ)
'RITE(‘.1050)
WRITE(‘.1060) (IPORTD(I.J).J=1.IPD(I))

400 CONTINUE

CLOSE(UNIT=5 . STATUS= 'IEEP ' )

RETURN
84D

fiffififif5

.DFIELDOOOOOOOOO......OOOOOO..........OOOOOOOOOOOOOOOOO00......

SUBROUTINE DFIELD (MD . MR. N1 .NP1 . NP2)

C“° DECLARATIONS
c

$INSERT syncnx
$INSERT CAUSBR
$INSERT UTILBK
slnsznr panmnx

c

C

C --- DFIELD GETS ALL INFORMATION FOR THE DYNAMIC ORSOURCE FIELDS
C AJDASHQT II) 1111-: R-FIELDS FROM THE ORIGINAL SYSTEM.

C

C INPUT --- MD. MR .N1.NP1.NP2 .

C OUTPUT --- MD.NR.N1.NP1.NP2 .

C (F THE MD JUNCTION as THE D-FIELD.

C

C L --- NODE INDEX

C 18 -- BOND INDEX

C IP --- PORT INDEX

C NUMJ --- BRANCH INDEX

C INDXDE --- NODE SERIEL NUMBER

C IDLSTD -- EBIENT TYPE LI ST

C IDNAM --- HEIENT NAME LIST

C NPTRD -- POINTER LIST FOR NBIMXD(START OF BOND GROUP)
C NBINXD -- LIST OF BONDS INCIDENT TO EACH EEMENT

C JUNCTD -- LIST OF BRANCHES IN THE D-FIELD

C LD --- NODE NUMBER STORAGE ARRAY FOR EACH D-FIELD
C IBD --- BOND NUMBER STORAGE ARRAY

C IPD -- PORT NUMBER STORAGE ARRAY

C

COOOCOOOOOOOOOOO0.3.0.0...OO‘OOOOOOOOOOOOO$00.0.........OOOOOOOOOOO

C
C --- INITIALIZATION
C

L=0

IB-IO

IP-O

40

45

C

106

NUMJ=0
GET JUNCTION DATA

L=L+1

NPTRD(MD.L)=IB+1

II=IB

DO 30 I=NP1.NP2

IF (ICMXT(I).EQ.0) RETURN
IB=II+I-NP1+1
NBIMXD(MD.IB)=NBIMX(I)
ICMXD(MD.IB)=ICMX(I)
IBOND=NBIMX(I)
IBMXD(MD.IBOND.1)-IBMX(IBOND.1)
IBMXD(MD.IBOND.2)=IBMX(IBOND.2)
CONTINUE

INDXDE(MD.L)=N1
IDLST(MD.L)=IELLST(N1)
IDNAM(MD.L)=IELNAM(N1)
IF(L.E0.1) GOTO 50

D0 40 I=1.L-1
IF(INDXDE(MD.I).EQ.INDXDE(MD.L)) GOTO 45
CONTINUE

GOTO 50

IB=IB+NP1-NP2

L=L+1

GOTO 85

GET ADJOINING NODE DATA

D0 80 I=NP1.NP2
IBOND=NBIMX(I)

NA= IBMX(IBOND.1)
NB= IBMX(IBOND.2)
N2=NA

IF(N1.EQ.NA) N2=NB

C-- FIND THE NODE WHICH HAS BEEN INCLUDED.

C

NUM=0
D0 55 J=1.20
IF(NBIMXD(MD.J).NE.IBOND) GOTO 55
NUM=NUM+1
IF(NUM.NE.2) GOFO 55
ICMXT(I)=O
GOTO 80
CONTINUE

D0 56 J-1.NBD2
IF (J.EQ.I) GOTO 56
IF(NBIMX(J).EQ.IBOND) GOTO 57

56

57

75

651
652

67

68

80
C

107

CONTINUE

IF(ICMXT(J).E0.0) GOTO 65
IF(IELLST(N2).GT.5) GOTO 59
ICMXT(I)=0

ICNXT(J)=0

L=L+1

IB=IB+1

NPTRD(MD.L)=IB
INDXDE(MD.L)3N2
NBIMXD(MD.IB)=NBIMX(I)
ICMXD(MD.IB)=ICMX(J)
IDLST(MD.L)= IELLST(N2)
IDNAM(MD.L)‘IELNAM(N2)
GOTO 80

D0 75 R=I.L
IF(INDXDE(MD.K).E0.N2) GOTO 80
CONTINUE

ICMXT(I)=0

NUMJ=NUMJ+1

JUNCED(NUMJ)=N2

JUNCB2(NUMJ)=I
JUNCB(NUMJ)'IBOND
JUNCOL(NUMJ)=N1

GOTO 80

ICNIT(I)=O

L=L+1

IB=IB+1

NPTRD(ID.L)=IB
INDXDE(MD.L)=N2
NBIMXD(ND.IB)-NBIMX(I)
IBOND=NBINX(I)

D0 651 181.NBD2

IF (1.80.1) GOTO 651

IF (NBIMX(J).EO.IBOND) GOTO 652
CONTINUE
ICNXD(MD.IB)=ICMX(J)
IPBIP+1

IPORTD(MD.IP)=IBOND
IF(IELLST(N2).EQ.6) GOTO 67
IF (IELLST(N2).EQ.7) GOTO 68
GOTO 80

IDLST(MD.L)=4
IDNAM(ND.L)='SE'

GOTO 80

IDLST(RD.L)=5
IDNAN(MD.L)='SF'

CONTINUE

C -- NEXT JUNCTION

C

108

NT=N1
85 IF(NUMJ.EQ.0) GOTO 120
N1=JUNCED(NUMJ)
NUMJ=NUMJ-1
NP1=NPTR(N1)
NP2=NPTR(N1+1)-1
DO 90 I=NP1.NP2
IF(ICMXT(I).EQ.1) GOTO 100
90 CONTINUE
GOTO 10

100 IXBR=IXBR+1
JBR(IXBR)=N1
N2=N1
L=L+1
IB=IB+1
NPTRD(MD.L)-IB
INDXDE(MD.L)=N2
NBIMXD(MD.IB)=NBIMX(JUNCBZ(NUMJ+1))
IBOND=JUNCB(NUMJ+1)
DO 653 J=1.NBD2
IF (J.E0.JUNCBZ(NUMJ+1)) GOTO 653
IF (NBINX(J).EQ.IBOND) GOTO 654
653 CONTINUE
654 ICIKD(MD.IB)=ICHI(J)
IP=IP+1
IPORTD(MD.IP)=NBINX(JUNCBZ(NUMJ+1))
IF(IELLST(N2).EQ.6) GOTO 110
IF(IELLST(N2).EQ.7) GOTO 115
GOTO 120
110 IDLST(MD.L)-4
IDNAM(MD.L)='SE'
GOTO 120
115 IDLST(MD.L)=5
IDNAM(MD.L)='SF'
120 IF(NUMJ.NE.0) GOTO 85
IPD(MD)=IP
LD(MD)=L
IBD(ID)=IB
NELD=L
IBD1=IB
NPTRD(MD.L+1)=IB+1
N1=NT

no

RETURN
END

C
C
C
COREARRG.0000......0000......CCOOOOOOOCCIOCCCCCOCOO00.000.00.00...t
C

SUBROUTINE REARRG(MRR.MDD.IB.NELS.IBMXS.NBIMXS.IBMXSN.

109

+ NBIXSN.INDXES)
C --- REARRG REARRANGE ALL THE ARRAY IN SEQUENCE.
C
C“‘ DECLARATIONS
C

INTEGER IA(60).IBB(60).NBIMXS(5.60).IBM(5.30.2).
+ IBMXS(5.30.2).NBIXSN(5.60).IBMXSN(5.30.2).INDXES(5.30)

c
$INSERT PARTBX
$INSERT syncs:
c
c
COOQOOOOOOOCOO00.0000.000.00.00.COCOOOOOOOCOOOOOCOOOO0.0.0.0000.0.0
c

IF(NRR.E0.0.AND.IDD.E0.0) RETURN

IF(MRR.EQ.0) GOTO 2

MS=MRR
XK=1
GOTO 5

2 MS=MDD
KK=2

C

5 D0 10 I=1 50
IA(I)=0
IBB(I)=0

10 CONTINUE

C

C
NBDSZ=IB
NBDS=NBDSZ/2

DO 100 I=1.NBD82
IA(I)=NBIMXS(MS.I)
100 CONTINUE

D0 120 I=1.NBDS2-1

LFI+1

D0 120 J=L.NBDSZ

IF(IA(J)-IA(I)) 110.120.120
110 ITEMP=IA(I)

IA(I)=IA(J)

IA(J)=ITEMP
120 CONTINUE

DO 130 I=1.NBDSZ

X=I/2.0

R=II2

YIX-K

IF(Y.EQ.0) IBB(X)=IA(I)
130 CONTINUE

DO 160 I=1.NBDS2
DO 150 J81.NBDS
IF(NBIMXS(MS.I).NE.IBB(J)) GOTO 150

150
160

170

180
190

193
194
192
191

200
210
220

110

NBIXSN(MS.I)=J
GOTO 160
CONTINUE
CONTINUE

DO 170 I-1.NBD
IBM(MS.I.1)'IBMXS(MS.1.1)
IBM(MS.I.2)=IBMXS(MS.I.2)
CONTINUE

DO 190 I=1.NBD

D0 180 J=1.NBDS
IF(IBM(MS.I.1).EQ.0) GOTO 190
IF(I.NE.IBB(J)) GOTO 180
IBMXSN(MS.J.1)=IBM(NS.1.1)
IBMXSN(MS.J.2)=IBM(MS.1.2)
CONTINUE

CONTINUE

DO 191 J=1.NBDS

D0 192 R=I.2

D0 193 I=I,NELS

IF (IBMXSN(MS.J.X).EQ.INDXES(MS.I)) GOTO 194
CONTINUE

IBMXSN(MS.J.K)=I

CONTINUE

CONTINUE

DO 220 N=1.5
IBOND=INTER(N.RX)

IF (IBOND.EQ.0) GOTO 220
R184

IF (RX.EO.1) [1'3

IF (INTER(N.K1).NE.MS) GOTO 220
D0 200 M81.NBDS
IF(IBB(M).EQ.IBOND) GOTO 210
CONTINUE

INTER(N.KX)=M

CONTINUE

RETURN
END

111

C
C—AUGIEN—OOOOOO0COCOOOOOOOOOOCOCOOOOOOOC00.00000COOOOOOOOOOOOOOOOOOOO
C

SUBROUTTNE AUGMEN(IR)

C -- PROGRAMMER -- TUNG ZHOU. AUG. 1983

C -- AUGMEN CONVERTS A STATIC SUBSYSTEM TO A DYNAMIC SUBSYSTEM BY
C INTRODUCING C OR I ELEMENTS FOLLOWING A GENERAL RULE.
C“‘DECLARATIONS

: DIMENSION NR(35).JZERO(10).JONE(10).IN10(35)

smsmr mam:
$INSERT PARTBX
smsmn‘ CAUSBK
$INSERT UTILBK
c

CO....OOOOOOOO.........O...‘0................OO.......OOOOOOOOOOOOCOO

C
C--- RESTORE.THE SUBSYSTEM DATA FOR PROCESSING.

C
NEL=LR(IR)
NBDZ'IBR(IR)
NBD=NBD2I2
C
CALL UNINAM(IR.IRLST.IRLNAM.NPTRR.IBMXRN.ICMXR.NBIXRN)
C
C-- ADD DYNAMIC ELEMENTS TO ACTING JUNCTIONS.
C

'RITE(.:5)

DO 10 J'I.INN
IF(INTER(J.3).NE.IR) GOTO 10
N2=IBMX(INTER(J.I).2)

IF (IELLST(N2).EQ.6.0E.IELLST(N2).EQ.7) JNOD=N2
IF(IELLST(N1).EQ.6.0R.IELLST(N1).EQ.7) JNOD=N1
IF (JNOLD.EQ.JNOD) GOTO 10
GOTO 15

C
5 FORMAT(/.'ADD C/I ELEMENT TO THE ACTING JUNCTIONS')
15 CALL ADELMT(JNOD.IR)

JNOLD=JNOD
C
10 CONTINUE
C
C
C---- TEST FOR COMPLETED CAUSALITY ON THE GRAPH.
C

400 D0 20 J=1.NBD2
IF (ICMX(J).EO.1) GOTO 30
20 CONTINUE

112

GOTO 800
C
C -- ADD C- OR I- ELEMENTS TO THE MULTI-R TYPE JUNCTIONS.
C

30 IRITE(‘.40)

40 FORMAT(/.'ADD C/I ELEMENT TO THE MULTI-R TYPE JUNCTIONS')
200 D0 50 131.NEL

50 NR(I)=0

DO 60 181.NEL
NP1=NPTR(I)
NP2=NPTR (1+1 )-1
IF (NP1.EQ.NP2) GOTO 60
DO 55 R=NP1.NP2
IBOND=NBIMX(K)
IF (ICMX(X).NE.1) GOTO 55
NAFIBMX(IBOND.1)
NB=IBMX(IBOND.2)
N2=NA
IF (I . EQ.NA) N2=NB
IF (IELLST(N2).EQ.3) NR(I)=NR(I)+1
55 CONTINUE
60 CONTINUE

JOLD=1
DO 70 J-1.NEL
IF (NR(J).LE.NR(JOLD)) GOTO 70
JOLD=J
70 CONTINUE

IF(NR(JOLD).EQ.1) GOTO 85
IF (NR(JOLD).E0.0) GOTO 800
JNOD=JOLD

CALL ADELMT(JNOD.IR)

GOTO 200

C
C
C -- ADD C- OR I- ELEMENTS TO THE.O- OR 1- JUNCTIONS WHICH
C POSESSES THE MOST INTERNAL BONDS.

C

85

'RITE(‘.90)
90 FORMAT(/.'ADD C/I ELEMENT TO THE MULTI-BRANCH JUNCTIONS')
95 DO 100 I=1.NEL

IN10(I)=0
100 CONTINUE

D0 120 N1=1.NEL

NP1=NPTR(N1)

NP2=NPTR(N1+1)-1

IF(NP1.EQ.NP2) GOTO 120

D0 105 X=NP1.NP2

IF (ICMX(X).EQ.1) GOTO 110
105 CONTINUE

110

115
120

130

320

330

340

345

310

113

OOTO 120

D0 115 K=NP1.NP2

IBOND=NBIMX(K)

N2=NA

IF(N1.EQ.NA) N2=NB

IF(IELLST(N2).EQ.6.0R.IELLST(N2).EQ.7) INIO(N1)=IN10(N1)+1
CONTINUE

CONTINUE

JOLD=0

DO 130 181.NEL

IF (IN10(I).LE.2) GOTO 130

IF (IN10(I).LE.IN10(JOLD)) GOTO 130
JOLD=I

CONTINUE

IF (JOLD.E0.0) GOTO 300
JNOD=JOLD

CALL ADELMT(JNOD.IR)
GOTO 95

ADD C- OR I- ELEMENTS TO THE SINGLEPR.TYPE JUNCTIONS.

J0=0

J1=0

IRITE(‘.301)

FORMAT(/.'ADD C/I ELEMENT TO THE SINGLEPR TYPE JUNCTINS')
D0 310 N=I.NEL

IF (IELLST(N).EQ.6.0R.IELLST(N).EO.7) GOTO 320
GOTO 310

NP1=NPTR(N)

NP2=NPTR(N+1)-1

D0 330 X=NP1.NP2

IF (ICMX(K).EQ.1) GOTO 340

CONTINUE

GOTO 310

IF(IELLST(N).EQ.7) GOTO 345

J0-J0+1

JZERO(J0)'N

GOTO 310

J1-J1+1

JONE(J1)=N

CONTINUE

J10'J1
IF(J1.GT.J0) J10=J0

D0 350 J31.J10
JONE(J)=JONE(J)
IF (J1.GT.J0) JWE(J)=JZERO(J)

114

NP1*NPTR(JONE(J))
NPZ-NP'I'R(JWE(J)+1)-1
D0 355 X=NP1.NP2
IF (ICMX(K).EQ.1) GOTO 360
355 CONTINUE
GOTO 350
360 JNOD=JONE(J)
CALL ADELMT(JNOD.IR)
350 CONTINUE
GOTO 800

800 D0 850 J=1.NBD2
IF (ICMX(J).EQ.1) GOTO 900
850 CONTINUE

CALL BARNAM(IR.IRLST.IRLNAM.NPTRR.IBMXRN.ICMXR.NBIXRN)

900 RETURN
END

C—UNINAM—OOOOOOO0.0.0.0.0....0.0.30...0........OOOOOOOOOOOOOOOOOO

C
SUBROUTINE UNINAM(IRD.ILST.INAM.NPT.IBM.ICM.NNBIM)

--- UNINAM STORES THE DATA OF A SUBSYSTEM INTE THE MATRICES
WHICH ARE COMPATIBLE WITH THEVPROGRAM ENPORT5(GLOBAL FORM).

(3C3f5()

DIMENSION ILST(5.30).NPT(5.31).IBM(5.30.2).
+ ICM(5.60).NNBIII(5.60)
cnmcmuz INAM(5.30)
smsmrr SYBGBK
INSERT mun:
stnsanr CAUSBK
gmsmr mus:

c...000......000............OOOOOOOOOOOO0.000000000000000000000...
.C

C

C
D0 10 I=1.NEL
IELLST(I)'ILST(IRD.I)
IELNAM(I)=INAM(IRD. I)

10 CONTINUE
DO 15 I‘1.NEL+1
NPTR(I)‘NPT(IRD.I)
15 CONTINUE

D0 20 I81.NBD
IBMX(I.1)=IBM(IRD.I.1)

115

IBMX(I.2)=IBM(IRD.I.2)
20 CONTINUE
C
DO 25 I=1,NBD‘2
ICMX(I)=ICM(IRD.I)
NBIMX(I)=NNBIM(IRD.I)
25 CONTINUE
C
RETURN
END
C
C
C
C—BAKNAu—OOOOC3000...0.00.0.0...00.0.0000000000000COCOOOOCOOOOOOCO
C
SUBROUTINE BAINAM (IRD.ILST.INAM.NPT.IBM.ICM.NBIM)
C
C --- BARNAM RESTORES THE DATA OF THE SYBSYTEM BACK TO ITS LOCAL
C MATICES.
C

DIMENSION ILST(5.30).NPT(5.31).IBM(5.30.2).ICM(5.60).
+ NBIM(5.30)
cnmcrmoaz INAM(5.30)
$INSERT SYBGBK
$INSERT PARTBX
$INSERT causnx
$INSERT mm:
c
c.0000.....OOOCCOOCOCOCCOOOOO‘COCOOOOOCOCOOOOOOCOC.30.000.00.000...0

C

C
C
NBD2=NBD‘2
DO 10 I=1.NEL
ILST(IRD.I)=IELLST(I)
10 INAM(IRD.I)'IELNAM(I)
C

DO 15 I=I.NEL+1
15 NPT(IRD.I)=NPTR(I)

C
no 20 I‘1.NBD
IBM(IRD.I.I)'IBMX(I.1)
20 IEM(IRD.I;2)3IEMX(Ip2)
C
ICM(IRD.I)'ICMX(I)
25 NEIM(IRD.I)‘NBIMX(I)
C
RETURN
END
C
C

C—mw-OOOO0.0...0.00.30...O...0.0.0.0..........OOOOOOOOOOOOOOOOO0

116

C
SUBROUTINE ADELMT(JNOD.IR)
C
C -- ADELMT ADDS THE.I OR C ELEMENT TO AN APPROPRIAT JUNCTION.
C

SINSERT SYBGBR
$$NSERT PARTBK
$INSERT CAUSBR

C
COOOOOOCOOOOOCOOCOCOOOCOO......OCCCOCO0.00.00.00.00000000QOOOCOOOOOO
C
C
IF(IELLST(JNOD).EO.6) GOTO 110
NENAMEs'I'
IBTD=1
GOTO 120
110 NENAME='C'
IBTD=-1
120 NEL=NEL+1
NBD=NBD+1

WRITE(‘.1000) NENAME.JNOD.INDIRE(IR.JNOD)
1000 FORMAT(/.' ADD ',A2.' ELEMENT‘TO NEW NODE '.I2,
+ ' (THE OLD NODE '.I2.' )')
IF(IBTD.EQ.-1) IELLST(NEL)=1
IF(IBTD.EO.1) IELLST(NEL)-2
IELNAM(NEL)=NENAME

IBMX(NBD.1)‘JNOD
IBMX(NBD.2)=NEL
N1=NPTR(JNOD)
N2=NPTR(NEL)-1
J8N2
DO 50 N=N1.N2
NBIMX(J+1)'NBIMX(J)
ICMX(J+1)=ICMX(J)
J=J-1

50 CWTINUE

NBIMX(N1)-NBD
NBIMX(N2+2)'NBD
NBIT=6
IF(IBTD.EQ.1) NBIT=3
ICMX(N1)'9-NBIT
ICMX(N2+2)=NBIT
NI-JNOD+1
DO 60 N8N1.NEL
NPTR(N)'NPTR(N)+1

60 CONTINUE
NPTR(NEL+1)'NPTR(NEL) +1

MNSTRN=40
MNSTRM=40

270
280

117

LEVEL=2
IF(IELLST(NEL).EQ.1) THEN
IT*6
ELSEIF (IELLST(NEL).EQ.2) THEN
IT=3
ELSE
GOTO 280
ENDIF
NP1=NPTR(NEL)
NP2=NPTR(NEL+1)-1

DO 270 X=NP1.NP2
ICMX(K)=1

CALL TWPASS(NEL.X.IT)
CONTINUE

CONTINUE

RETURN
END