Michael Balchanos, Santiago Balestrini, Dr. Neil Weston, and Dr. Dimitri Mavris

Multi-Physics Time-Variant First-Order Model Integration of Complex Systems ABSTRACT This paper introduces one of the key areas of research being conducted at the Aerospace Systems Design Laboratory (ASDL) in support of the Integrated Engineering Plant (IEP) initiative (ASDL 2004) funded by the Office of Naval Research. The intent of this effort is to develop a method to integrate heterogeneous multi-physics dynamic models into a single framework to produce a comprehensive modeling and simulation environment of complex systems. The key challenge consists of accurately cosolving time-domain simulations in an efficient manner. Such systems, similar to the Integrated Engineering Plant (IEP), typically include systems that fall within multiple physical categories (Dunnington, Garter, and Stevens 2003), e.g., complex systems that incorporate mechanical, hydraulic, electrical or other subsystems. The need for a method to facilitate software model integration is briefly introduced, followed by a summary of previous efforts to integrate heterogeneous dynamic systems. The background section contains information on supporting techniques that will prove to be enablers for the methodology presented here. The methodology is introduced and described in detail, in order to enable readers to apply it to their application. Finally a proof-of-concept is presented with the accompanying results to demonstrate the feasibility of employing this methodology to real world problems. The results are divided into two sections. The proves that the previous method, referred to as the zeroorder integration method, is proven to generate results that are counter to what would be expected and could mislead simulation integrators. The second portion analyzes the difficulties in implementing the first-order

method and provides insight as to where these arise from.

INTRODUCTION US Navy ships are becoming more complex as more systems are integrated in order to increase mission effectiveness and reduce ownership costs. As the complexity of these systems increases, the ability of monitoring and controlling their behavior is reduced. Main reasons for this difficulty are the nonlinear behaviors (disproportionate cause-to-effects), the multitude of interactions and the heterogeneity of the various subsystems, which complicate the understanding of the behavior of the system and make predictions on their behavior a challenging task. Due to their complexity, understanding their behavior can only be done through modeling and simulation. A major complication in the process of simulating integrated systems is that they span multiple realms of science and engineering (electromagnetism, fluid mechanics, thermodynamics, controls, etc) (Hughes et al. 2006). The need to model this multitude of disciplines as accurately as possible forces the need to bring together knowledge from different disciplinarians, in order to obtain a comprehensive model that captures the integrated behavior of the system. To accomplish this task there are two main routes that can be pursued. One route is to use models that are created by domain experts and can be programmed into a single modeling and simulation framework. The second route is to allow experts to model the various subsystems using the framework of their choice and then integrate the models a posteriori. The first option is generally denominated direct translation (Monti 2006) because the models are translated to the single modeling framework of choice. It offers some benefits over the second option in that it produces a model that is better

© 2007 by Balchanos, Balestrini, Weston, and Mavris. Published by the American Society of Naval Engineers with permission.

behaved (time domain integration issues are avoided) and it tends to run faster.




discipline oriented software model developers (disciplinarians) have preferred modeling languages or software packages for model development,




at times validated and accepted models already exist in a given framework and there is no easy access to the mathematical equations behind the model,




not using the tools preferred by domain experts can hinder acceptance of the integrated model, and




as models are updated by the domain experts, the entire integrated model needs to be updated.

Thus, many times it may be desirable to a posteriori integrate time-domain models of complex systems from multi-disciplinary fields. It is important to note that these models have independent time solvers, and cannot be solved by a single high-level solver. Currently, there is a gap in the capability to accomplish this. This paper will discuss a process for accomplishing this based on existing methods that have been developed to intelligently integrate dynamic software/hardware models, in conjunction with the use of Bond Graph Theory to generalize the method.

Previous Work A zero-order integration approach was previously developed, using a technique that has been introduced as the Multi-Disciplinary Simulation (MDS) (Weston et al. 2006). This method proved to be useful (Hughes et al. 2006) but highlighted several problems within timedomain integration. The problem results from the fact that fidelity (how well the analysis represents reality) has to be traded off with demand for computational effort, alternatively referred to as the fidelity-efficiency tradeoff. Figure 1 depicts how by applying a new technology, the tradeoff between effort and fidelity is expected to be closer to the ideal.

Zeroth Order Method

Fidelity

Nevertheless, it has its drawbacks; the main ones being that

First Order Method Goal

Computational Effort

FIGURE 1. Tradeoff between fidelity and effort

Experience with a preliminary IEP model revealed that it was necessary to integrate all subsystem models using the smallest time step that was required by the subsystem with the fastest dynamics, in order to maintain accuracy and stability. Furthermore, in some instances, the simulation yielded completely inaccurate results that greatly misrepresented reality. This motivated the research team to pursue an approach that attempts to capture the interactions between the models more accurately to reduce this error, while at same time allowing each subsystem model to be integrated with the most appropriate time step.

BACKGROUND Energy-based software integration Method The application of energy conservation methods to develop the governing equations associating the states of different physical subsystems is a common practice in science and engineering. Despite the fact that multidisciplinary systems consist of heterogeneous physical systems, they all exchange energy amongst each other, of the rate of exchange of this energy is power. In other words, power plays the role of the “global” variable that can be estimated for any physical system. Energy conservation will ensure that power that is exchanged between two heterogeneous systems (e.g., electrical with thermal system) is the same, yet the variables that are used to describe this power for every system are usually different. Even though the concept of conservation of power spans all fields of physics, it is necessary to identify a common notation for all. After an extensive literature review, Bond Graph Theory (BGT) was selected as an appropriate enabler for this purpose because it provides guidelines for

© 2007 by Balchanos, Balestrini, Weston, and Mavris. Published by the American Society of Naval Engineers with permission.

abstracting and modeling the dynamics of any combination of heterogeneous physical subsystems.

e0 = F0

Bond Graph Theory (BGT)

f0 = u0

BGT is a tool used to analyze and describe dynamical systems and provide a generic graphical process that allows for the derivation of the governing equations for any combination of any type of dynamical system. BGT relies on the concept of the tetrahedron of state, depicted in Figure 2, and recognizes that the relationship between these states is independent of the physical domain on which it is applied. The theory is based on the universality exhibited by all dynamical systems. Figure 3 provides a representation of this phenomenon, along with the Bond Graph representation of three systems that are dynamically identical. e

e1 = F1 m1 f1 = u1

b1 (a) Mechanical System

e0 = v0

m1

b1

f 0 = i0

f 1 = i1

k1

(b) Electrical System

effort flow momentum displacement

∫d

t

C

e: f : p: q:

k1

f 1 = Q1

R

f0 = Q0 q

e1 = P 1

e0 = P0 (c) Hydraulic System

d/d

t

p I

f

R : resistance I : inertance C : capacitance

FIGURE 2. The tetrahedron of state

BGT is an enabler for modeling and describing complex dynamical systems. A Bond Graph is an explicit graphical tool for capturing the common energy structure of physical systems (Brown 2001). Such physical systems are represented as interconnected symbols with lines, denoting power flow paths. The energy structure is represented by a generalized power P that is expressed as P =e⋅ f

(1)

where e is an effort and f is a flow rate. The two latter variables are generalized representations of more problem specific variables. For instance, in

e0 f0

C:k1

I:m1

e0 fk

em f1

0

e0 f1

1

e1 f1

eb f1

R:b1 (d) Bond Graph Representation FIGURE 3. The universality of dynamical systems, adapted from Gawthrop and Bevan (2007)

© 2007 by Balchanos, Balestrini, Weston, and Mavris. Published by the American Society of Naval Engineers with permission.

the case of an electrical system where voltage v and current i are state variables of interest, the corresponding Bond Graph would include electric power P as the product of voltage times the current. The voltage v corresponds to an effort and the current i corresponds to a flow. Table 1 lists the effort and flow of each physical domain. TABLE 1 Efforts and flows by discipline Physical Domain

Effort ( e )

Flow ( f )

Electrical

Voltage (v)

Current (i)

Hydraulic

Pressure (p)

Thermo-Fluid

Temperature (T)

Mechanical

Force (F)

Mechanical Rotational

Torque (T)

Volume Flow Rate (dV/dt) Entropy Flow Rate (ds/dt) Velocity (v) Angular Velocity (ω) Mole Flow Rate (dN/dt)

Chemical

Chemical Potential (µ)

Thermo-Chemical

Enthalpy (h)

Mass Flow Rate (dm/dt)

Magneto-motive Force (em)

Magnetic Flux (φ)

Magnetic

Complex physical systems are interconnected by bonds/junctions to form a flow network structure. Subsequently, the system equations can be systematically derived by the power flow network structure. A typical example of the aforementioned process is shown in Figure 4, for the purpose of deriving the governing differential equation for a spring-mass-damper system. Free-body Diagram

Bond Graph Representation

R:b uref = 0 k

b

1/k .. C

uref 1

I:m

SE : mg m uref = 0 mg

Differential Eq. Representation

d 2x dx m 2 + b + kx = − mg dt dt

FIGURE 4. Application of Bond Graph Theory to a spring-mass-damper system, based on figure in p. 86 Karnopp, Margolis and Rosenberg (2006)

Time-Variant First-Order Approximation Method A zero-order interfacing scheme, dubbed as the Multi-Disciplinary Simulation (MDS), was developed and implemented (Weston et al. 2006) for the purpose of software model integration in order to produce dynamic simulations of large complex systems. However, several issues have appeared, in particular the poor tradeoff capability between fidelity and effort as depicted in Figure 1. The goal is to reduce computational effort by extending the simulation time steps of subsystems with larger time constants and vice versa. An alternative research direction that has been investigated for improving the integration of heterogeneous software models was inspired from a method for interfacing Power-Hardwarein-the-Loop (PHIL) simulations for simple R-L or R-C electric circuits (Wu 2005). PHIL simulations are an extension of Hardware-in-theLoop (HIL) techniques, which additionally enable the natural coupling (continuous state variables, e.g., voltage, current, etc.) in addition to the signal coupling (discrete variables, e.g., on/off, etc.) between the simulation environment and the system hardware. In other words, this interfacing scheme includes energy conservation, so that real power is virtually exchanged between simulated and physical power components. According to this scheme, the integration of the Hardware Under Test (HUT) with the Rest Of the System (ROS) is expressed as a first order differential equation. The form of this equation can generally be for an inductive/resistive combination: di HUT = α ⋅ i HUT + β ⋅ v ROS dt

(2)

or for a capacitive/resistive combination: dv HUT = α ⋅ i ROS + β ⋅ v HUT dt

(3)

The extension of this method to the integration of software-software simulators has been identified as a possible solution to the challenge presented in this paper. Nevertheless, the assumption of capacitive/resistive and

© 2007 by Balchanos, Balestrini, Weston, and Mavris. Published by the American Society of Naval Engineers with permission.

inductive/resistive may not be applicable to some of the subsystems of interest. Furthermore, it is critical to analyze the stability of the integrated system solution to the extrapolation of rates of change of interface states. This means that depending on how far the states are extrapolated, i.e., the time step used to coordinate the subsystem simulations, the numerical solution of the integrated system may exhibit divergent behavior to a greater or lesser degree.

DEVELOPMENT OF THE MODEL INTEGRATION METHODOLOGY The ultimate purpose of this methodology is to provide guidelines that will allow the integration engineer to construct an integrated software model of a complex system. It is assumed that the software simulations of the system components are provided. The essence of this effort is concentrated on interfacing the various models of the system components in such a way that the behavior of the integrated system represents the true behavior of the real integrated system as precisely as possible.

Step 1: Understanding the Systems The purpose of this step is to define the complex system that will be created by the software integration from the two or more dynamical system simulations. The initial effort concentrates on understanding the interfaces of these, what they are, and how they are related. This method does not assume that it is necessary to know the internal operations of the models being integrated since in many cases this is not possible. Nevertheless, a general concept of how the models behave is important to enable the integrator to tailor the process to the specifics of the models of interest. System engineering tools such as the Tree diagram or the Affinity Diagram can be used to carry out the physical and functional decomposition of the complex system if its size is sufficiently large. For relatively small systems, simple diagrams can suffice. The result of this first exercise is a clear image of the interfaces of the component models that

will be integrated. These interactions are represented by the variables that describe the interactions and the assumptions inherent in the models to be integrated. As a final note, it is imperative that the explicit assumptions from the models be documented and understood as best as possible.

Step 2: Standardizing the Interactions Heterogeneous complex systems consist of subsystems that operate in disparate physical domains, e.g., mechanical, electrical, hydraulic, thermal, etc. For instance one model may represent a full scale electrical system and another might be a network of chillers, pipes and pumps to provide cooling for the electrical components. Nevertheless, it is the same energy that flows through two or more heterogeneous subsystems; however there are different variables in each system interface that are physically meaningful and are therefore used to evaluate the flow of energy. BGT was introduced in the previous section as a generalized framework for analyzing dynamical systems. The underlying philosophy of BGT, that power is conserved through multi-physics systems through the exchange of efforts and flows, is useful when determining how the different models are to interact. By using the standardized definition of efforts and flows, integrators can establish a clear dialogue with the domain experts developing the models. For a list of what these efforts and flows are, refer to Table 1. The result of this step is a standard representation (Bond Graph efforts and flows) of the interfaces between the models.

Step 3: Integrating the Models This step centers on the explicit integration of the subsystem simulations. For the integration to be effective, the simulations being integrated must be able to be paused at predefined times, its states queried and updated, and the simulation restarted. Otherwise, the task of integrating these simulations is considerably more difficult. The simulations are then integrated in a framework that will allow for the automatic execution and the exchange and storage of data. This can be done in a commercial integration

© 2007 by Balchanos, Balestrini, Weston, and Mavris. Published by the American Society of Naval Engineers with permission.

framework (e.g., Phoenix Integration’s ModelCenter®, Engineous’ iSight, TechnoSoft’s Tool Integration Environment, etc.), a scripting language (e.g., Python, VBScript, C#, etc.), or a programming language (e.g., Java, C++, etc.).

been developed. Figure 5 contains the physical depiction of the system and the monolithic implementation in Simulink®. In

Seawater

VDC

Two critical considerations are at the core of the third step. One is the scheduling of the simulation, i.e., determining when each individual simulation will run. The simulations can be made to run for a fixed time step, variable time step, or differential time step. The first means that all sub-system simulations run for the same period of time and exchange data at this predefined times. For the second, the time step is allowed to vary depending on how the interface variables vary. Therefore as the system approaches steady state, the time step increases. The third option allows different simulations to run at different time steps. The second and third options lead to the most significant time savings and can be improved by the first-order method. By estimating the first-order behavior of the simulations with the larger time constants, the simulations with the shorter time constants can receive more accurate estimates of the interface states, leading to more accurate overall behavior. The second consideration concerns the synchronization of the interfacing states. It is crucial that the simulations read the states from the appropriate time step. For this reason, the data should be centralized in the integration framework and the subsystem simulations should not be allowed to exchange states directly, but only through mediation of a synchronization repository that will ensure that the states read match their corresponding time.

DEMONSTRATION OF THE METHODOLOGY The starting point for the implementation of the integration method is the understanding of the physical models that have to be combined together in order to create a complex system simulation model. The individual models are treated as “black boxes” in the sense that they perform calculations using their own solvers and integration time steps. For this specific demonstration, a dynamic model of a small scale (miniature) Integrated Engineering Plant has

Load

HX

Cold Plate (T)

VDC

DC Bus

PCM Pump VDC

PGM Out

VAC

Seawater

(a) Physical Representation

(b) Simulink® Representation FIGURE 5. Integrated Thermoelectric System

The transient model of the miniature IEP has been divided into three subsystems, namely an electrical power system, a thermal system, and an electro-hydraulic. Figure 6 is the graphical depiction of the miniature IEP model as a set of three interfacing “black boxes”. Note that the interfaces adhere to the efforts and flows described by BGT. The electrical system consists of an AC power generator that includes first order exciter dynamics (control of voltage for the operation of the AC power generator) along with a governor for adjusting the AC frequency. An AC/DC converter is used to produce DC voltage to support the operation of the service load and the

© 2007 by Balchanos, Balestrini, Weston, and Mavris. Published by the American Society of Naval Engineers with permission.

two pumps that regulate the flow for the cooling systems of the electro-hydraulic subsystem.

Electro-Hydraulic (R Load, Pumps) of

y e p tur tr o era En of mp e Te rat ge an low Ch eF lu m

te Ra

Cu rre nt

Vo lt

ag e

Vo

re ssu Pre

Electrical (PGM, PCM, Bus)

Rate of Change of Entropy Temperature

testing the proposed integration method. Figure 7 is a screen capture of the model integrated in ModelCenter®. The Synch script is the data repository that ensures that the subsystem models query the appropriate values as described in Step 3 of the methodology. The scheduling of the simulation is implemented with the scheduler feature of ModelCenter®. R : FW

f1 : IFW

VDC : e0

Thermal System (Chiller, Heat Exchanger, Pipes)

e0 : VDC f0 : IDC

FIGURE 6. Interfaces of the integrated model

The thermal model consists of a closed-loop cooling system for extracting heat from the heat exchanger and an open-loop sea water system to cool down the closed-loop system through a chiller heat exchanger. The closed-loop cooling system includes a cold plate as the interface between the load and the heat exchanger as shown in Figure 5. A proportional thermostat is used for controlling the operation of the heat exchanger based on the temperature of the load. Heat exchange between the cold plate and the surrounding air is modeled. Figure 5 depicts this system.

e0 : VDC

0

SE

VDC : e0

R : Load

f2 : ILOAD f3 : ISW

R : SW (a) Electrical Interface.

R : Electrical

f1 : ds/dtElectrical

THX : e0 e0 : THX SE f0 : ds/dtTotal

0

e0 : THX

R : Load

f2 : ds/dtLoad

(b) Thermal Interface.

R : Pump FIGURE 7. Integrated system in ModelCenter®

The purpose of the monolithic model is to be the reference to which the integrated model will be compared. ModelCenter® by Phoenix Integration has been selected as an appropriate integration environment. The responses of the integrated environment will be compared to the corresponding responses of the monolithic (all in Simulink®) validated version of the model for

∆PPump : e1 e0 : ∆P SF f0 : dV/dtFW

f0 : dV/dtFW

1

e2 : ∆PPiping f0 : dV/dtFW

(c) Hydraulic Interface. FIGURE 8. Bond graph of interfaces

© 2007 by Balchanos, Balestrini, Weston, and Mavris. Published by the American Society of Naval Engineers with permission.

R : Pipes

Following the procedure for testing the integration approach, flow and effort variables have been identified for the employed models. Figure 8 displays the Bond Graph of the interfaces of the three subsystem models previously presented in Figure 6. These are a means to ensure that generalized power is conserved through the interfaces. A first-order integration method was developed, based on the TFA interfacing method. While in zero-order integration only a state variable (only efforts or flows, e.g., only voltage v or current i) would constitute the link between two different simulations, in the TFA implementation, there is also a rate (e.g., rate of the temperature T or rate of change of entropy ds/dt) that is being exchanged between models. The calculation of the rates of the flows and efforts at the interfaces for a given time step is being performed according to equations (2) and (3).

TABLE 2 Normalized Integration Metrics

310 Monolithic Integrated (TS 0.25 sec) Integrated (TS 0.50 sec) Integrated (TS 2.00 sec)

Temperature of Heat Exchanger (K)

305 300

Generalized Power

295 290 285 280 275

Entropy Flow Rate from Load Model (J/s-K)

270

0

50

100

150

200

250 300 Time (s)

350

400

450

500

2

1.5

1

0.5

implemented as a zero-order integration scheme. According to this method, only zero-order information is being propagated through the models, e.g., only state variable values, not their rates of change (first time derivative of the state variables), or the second time derivatives. Figure 9 demonstrates that the zero-order method is not suitable for integrating this dynamic model. Not only is the integrated behavior inaccurate, but the relation between fidelity and computational effort run counter to what intuition would dictate. Note, that the larger the time step, the better the behavior of the integrated system, i.e., the shorter the time step, the worse the integrated system does. This is a significant realization that needs further studying. Preliminary investigations indicate that as the time step increases, the interfacing through efforts and flows may actually lead to nonconserved power flow.

0

50

100

150

200

250 300 Time (s)

350

400

450

FIGURE 1. Comparison of THX for zero-order method with different time steps

RESULTS It was mentioned previously that the MDS method had exhibited several issues as it was

500

Run Time (s) Per Time Step Overhead GA Electric GB 2 (×10 ) GC GA Thermal (Electric) GB (×103) GC GA Thermal (Load/Pumps) GB (×10) GC GA Hydraulic (Fresh Water) GB (×10) GC GA Hydraulic (Sea Water) GB (×10) GC

Time step 0.25 sec 14193 7.09 2739% 5.71% 1.08% 6.97% 0.18% 0.18% 0.26% -6.81% 1.13% -0.36% -2.70% 1.72% -1.26% -42.55% -5.97%

Time step 0.50 sec 6944 6.94 1289% 4.41% 0.59% 6.03% 0.13% 0.13% 0.23% -4.87% -1.54% 0.57% -0.86% -11.83% 2.36% -3.20%

Time step 2.00 sec 1915 7.63 283% 0.15% -0.35% 5.07% -

The final analysis of results consists in the evaluation of the integration approach. The metrics defined in equations 4, 5, and 6 have been considered for the evaluation of the method. Measures f(t) and g(t) represent power responses from the monolithic model (prototype) and the integrated model in the ModelCenter environment (tested integration method) respectively. Evaluation metric values are summarized in Table 2 for the five generalized powers interfaced in the model.

© 2007 by Balchanos, Balestrini, Weston, and Mavris. Published by the American Society of Naval Engineers with permission.

t final

t initial

GB =

∫

f (t ) − g (t ) dt (4)

t final − tinitial

t final

t initial

f (t ) − g (t ) dt (5)

t final − tinitial t final

GC =

∫ [ f (t ) − g (t )]⋅ t dt (t − t ) t initial

(6)

2

final

initial

The metrics calculated using equations 4, 5 and 6 attempt to capture different errors produced by the integrated model. Metric GA, equally penalizes divergence, but it can cancel the error if the values alternate between positive and negative. The second metric, GB, does not cancel out if the error values alternate between positive and negative, but equally penalize if the error is at the beginning of the time history or at the end. The third metric, GC, weighs the steady state error more heavily than the initial transient error that generally occurs in integrated systems. The (-) in Table 2 represents zero and is the minimum error for that metric and generalized power, the other numbers of each row specify how much additional error that integration had with respect to the minimum error. As shown in Table 2, the larger time step integration presented the least amount of error for 12 of the 15 metrics computed.

Hydraulic Power of Fresh Water Loop (W)

0

-20

-40

-60

-80 Monolithic Integrated (TS 0.25 sec) Integrated (TS 0.50 sec) Integrated (TS 2.00 sec)

-100

-120

0

50

100

150

200

250 300 Time (s)

350

400

450

500

FIGURE 10. Generalized Hydraulic Power for Fresh Water Loop

The runtimes for the different integrations were computed and are presented at the top of Table

2. The overhead was computed as the difference between simulation time and runtime normalized to simulation time. In real-time, overhead would be zero and at half the speed of real-time it would be 100%. It is important to realize that the runtime per time step remains almost constant for the three integrations, even though in the subsystem solvers, the computations are for very different time spans. This exemplifies how much of the actual run time is devoted to simulation overhead (i.e., initialization, data swapping, etc.). The results in Table 2 are counterintuitive since they contradict the general shape of the notional Pareto front presented in Figure 1. The inverse of the error values presented in Table 2 are plotted against the total computational time in Figure 11 to support the previous observation. 1

Normalized Fidelity (Inverse of Error)

GA

∫ =

0.95

Metric A Metric B Metric C

0.9 0

5000

10000

15000

Computational time (sec)

FIGURE 11. Fidelity vs. Computational Effort for zero-order method

Of particular interest is the coupling of this result with the runtimes for each simulation, which would indicate that more fidelity can be attained with one tenth of the overhead. This result is not only counterintuitive but proves that integrating these time-domain simulations can be a daunting task. The simulation integrators do not have the luxury of having a monolithic model to compare their results with, at least not for every design. Validation of the integrated simulation environments cannot be done fully, and at best can only be done for a few settings of the design parameters or by parts. What this result illustrates is that validation by parts would prove to be erroneous, since even though the parts are the same for both systems, the integrated behavior is very different. Furthermore, if an integrator where to test

© 2007 by Balchanos, Balestrini, Weston, and Mavris. Published by the American Society of Naval Engineers with permission.

Regardless of the analysis of these results, it is clear that the zero-order method is not suitable for this integration. The authors estimated that a first-order method may be more successful. For this reason the TFA method was attempted for an interface that was considered critical to matching the monolithic and integrated simulations. Capacitive System

Rate of Change of Pressure Drop (Pa/s)

15

10

10

10

5

10

the method. Figure 11 presents the actual rate of change of the effort and flow in black and the estimated one by TFA in red. The behavior diverged considerably as depicted in Figure 12; note the scale of the rates is logarithmic. The root of the problem was traced to the determinant that must be calculated to invert the matrix. When this rate becomes zero, the matrix becomes singular, and can therefore not be inverted. Furthermore, the sensitivity of the system to these values, and the errors in the estimation of the gradients by TFA is clear evidence that the propagation of error would lead to considerable integration errors. 400

Rate of change of Pressure Drop (Pa/s)

different time steps to see what the behavior of the integrated system truly was, the results would point in the opposite direction. This would be because if it were assumed that a smaller time-step would be more accurate, and the trend extrapolated, the integrator would be forced to pick a result that would point to the opposite direction.

200 0 -200 -400 -600 -800 -1000 0

0

10

-2000 -4000 Pressure Drop (Pa)

0

50

100

150

200

250 300 Time (s)

350

400

450

500

-6000

0

0.002

0.006

0.004

0.01

0.008

0.012

0.014

Volume Flowrate (m3/s)

-4

x 10

Inertial System

20

10

15 Rate of change of Flowrate (m3/s 2)

Rate of Change of Volume Flowrate (m3/s 2)

10

5

10

0

10

-5

10

10

5

0

-5 0.015

-10

10

0

50

100

150

200

250 300 Time (s)

350

400

450

500

0.01 0.005

FIGURE 12. Rates of change of effort and flow

The variations in pressure and flow rate were considered critical due to their reciprocal sensitivity and to the overall solution. The pressure drop and the volume flow rate of the fresh water loop were chosen as the effort and flow respectively for which a first-order method based on Wu’s TFA would be implemented. The initial results of the integration showed that the first-order approximation did not improve

Volume Flowrate (m3/s) 0

0

-1000

-2000

-3000

-4000

-5000

-6000

Pressure Drop (Pa)

FIGURE 13. Trajectories of the rate of change of pressure drop and volume flow rate in fresh water loop as a function of pressure drop and flow rate

More extensive analysis of the instability points showed that the rates of change as a function of time are zero in these indeterminate zones, but the rates of change as a function of the state variables becomes infinite. This is the relation that concerns TFA since the rates are solely

© 2007 by Balchanos, Balestrini, Weston, and Mavris. Published by the American Society of Naval Engineers with permission.

calculated as a function of the effort and flow state variables. The red dots in Figure 13 correspond to the times that the TFA matrix becomes nearly singular.

CONCLUSIONS A first-order integration method has been developed and tested for integrating heterogeneous software models for the purpose of performing time-domain simulations of complex systems. Results comparing a monolithic model representing a miniature IEP system with an integrated version of the same model are presented and analyzed to understand the lost fidelity due to the integration of independent simulations. Evaluation of the analysis results has revealed a valuable set of lessons to consider when attempting to integrate simulations with independent solvers. Despite the fact that the approach seems to be headed in the right direction from a theoretical standpoint, in practice the convergence of the results is not as of yet satisfactory. The TFA method, in the form that it has been applied for this purpose, is not adequate for guaranteeing a robust integration environment. Further research must be conducted to address the discrepancies between the behavior of the monolithic and the integrated systems. In particular the lags observed. Additionally, the method requires the solution of a 2x2 system, but the presence of singularity points in this solution is responsible for an intensive unsteady behavior in the responses of the integrated system. Bond Graph Theory, originally identified as a suitable enabler for an energy-based integration method, seemed to prove to be inadequate when interface variables lagged. Since models are interfaced by a balance of effort and flows under the precept of conservation of power, the lag produces an inconsistency in the conservation of power, ultimately leading the integrated model to diverge from its monolithic implementation.

under contract N00014-07-1-0426. We want to further extend a very special thank you to Mr. Anthony J. Seman and Mr. Frank T. Ferrese for their invaluable suggestions and insightful comments which proved to be a great help in guiding this effort. Further appreciation is extended to Mr. George Stefopoulos for his generous assistance in the development of the high fidelity electrical model.

REFERENCES Aerospace Systems Design Laboratory (ASDL), “Integrated Reconfigurable Intelligent Systems,” Georgia Institute of Technology, 2004. Brown, F.T., “Engineering System Dynamics: Control Engineering,” 1st ed., Marcel Dekker, Inc., New York, NY, 2001. Dunnington, L., Garter, G. and H. Stevens, “Integrated Engineering Plant for Future Naval Combatants Technology Assessment and Demonstration Roadmap,” Systems Engineering Group, Engineering Technology Center, Marine Technology Division, Anteon Corporation. 2003. Gawthrop, P.J., and G.P. Bevan, “Bond-Graph Modeling: A Tutorial Introduction for Control Engineers,” IEEE Controls Systems Magazine, vol. 27(2), pp. 24-45, April 2007. Hughes, R.; Balestrini, S.; Kelly, K.; Weston, N.R., and D.N. Mavris, “Modeling of an Integrated Reconfigurable Intelligent System (IRIS) for Ship Design,” Ships & Ship Systems (S3) Technology Symposium Change, Challenges & Constants, West Bethesda, MD, November 13-14, 2006. Karnopp, D.C.; Margolis, D.L., and R.C. Rosenberg, “System Dynamics: Modeling and Simulation of Mechatronic Systems,” 4th ed., John Wiley & Sons, Inc., Hoboken, NJ, 2006. King, R.P., “Introduction to Practical Fluid Flow,” Elsevier Science, Oxford, 2002. Kundur, P., “Power System Stability and Control,” McGraw-Hill, Inc., NY, 1994. Monti, A., “Model Wrappers,” presented at the 2006 VTB Users and Developers Meeting, University of South Carolina, Columbus, SC. Retrieved May 31, 2007 from Mosterman, P.J., “Hybrid Dynamic Systems: A hybrid bond graph modeling paradigm and its application in diagnosis,” Ph.D. Thesis, Vanderbilt University, 1997. Sauer, P.W., and M.A. Pai, “Power System Dynamics and Stability,” Prentice Hall, NJ, 1998.

ACKNOWLEDGMENTS

Zivi, E.L., and R. Youngs, “Fluid-cooled power component heat exchanger model,” whitepaper, June 2005.

The authors wish to acknowledge the Office of Naval Research for sponsoring this research

Vasquez, H., and J.K. Parker, “A new simplified mathematical model for a switched reluctance motor in a

© 2007 by Balchanos, Balestrini, Weston, and Mavris. Published by the American Society of Naval Engineers with permission.

variable speed pumping application,” Mechatronics, vol. 14(9), pp. 1055-1068, November 2004. Weston, N.R., Balchanos, M.G., Koepp, M.R., and D.N. Mavris, “Strategies for Integrating Models of Interdependent Subsystems of Complex System-of-Systems Products,” IEEE Proceeding of the Thirty-Eighth Southeastern Symposium on Systems Theory, pp. 310-314, March 5, 2006. Wu, X., and A. Monti, “Methods for Partitioning the System and Performance Evaluation in Power-Hardwarein-the-Loop Simulations – Part I,” proceedings of the 31st Annual Conference of IEEE, pp.251-256, November 6-10, 2005. Wu, X., “Methods for Partitioning the System and Performance Evaluation in Power-Hardware-in-the-Loop Simulations,” PhD Thesis, University of South Carolina, 2005. Yoo, K.; Simpson, K.; Bell, M., and S. Majkowski, “An Engine Coolant Temperature Model and Application for Cooling System Diagnosis,” SAE 2000 World Congress, Detroit, MI, March 2000.

Mr. Michael Balchanos obtained his Diploma in Physics from the Aristotle University of Thessaloniki, Greece and his M.S. in Aerospace Engineering degree from Georgia Tech. Mr. Balchanos’ Ph.D. dissertation topic will entail the development of a robust methodology for the design of complex naval systems under mission effectiveness requirements, utilizing physicsbased modeling and simulation environments. He has also been involved on understanding and developing electrical models for simulating naval ship operations and contributing in the integration of time-domain simulations in a comprehensive integrated framework.

experience working on the different phases of design, analysis, and testing of various complex systems. His main areas of interest include developing and implementing new methods for systems design and optimization, developing physics-based simulation models, and all areas of theoretical, applied and computational mechanics. Dr. Dimitri N. Mavris is the Boeing Professor of Advanced Aerospace Systems Analysis at the Guggenheim School of Aerospace Engineering, Georgia Institute of Technology, and the director of its Aerospace Systems Design Laboratory (ASDL). Since his initial appointment as academic faculty in 1996, Dr. Mavris has consistently acted to expand his academic and teaching abilities and responsibilities. His primary areas of research interest include: advanced design methods, aircraft conceptual and preliminary design, and air-breathing propulsion design.

Mr. Santiago Balestrini obtained his B.S. and M.S. in Aerospace Engineering degrees from Georgia Tech. Mr. Balestrini’s Ph.D. dissertation topic is to develop a modeling process to understand complex systems. Under the IRIS initiative, Mr. Balestrini has helped integrate time-domain simulations in a comprehensive integrated framework that will help prototype hardware as well as control algorithms for the US Navy, effectively reducing cost and improving mission effectiveness. Dr. Neil R. Weston is a Research Engineer II with the Aerospace Systems Design Laboratory (ASDL) at the School of Aerospace Engineering of the Georgia Institute of Technology. Dr. Weston has more than fifteen years of © 2007 by Balchanos, Balestrini, Weston, and Mavris. Published by the American Society of Naval Engineers with permission.