Olivier Dalle
email: olivier.dalle@sophia.inria.fr

OSA : AN OPEN COMPONENT-BASED ARCHITECTURE FOR DISCRETE-EVENT SIMULATION

This paper describes work in progress to initiate the collaborative development of a new software platform for discrete-event simulation studies, the Open Simulation Architecture (OSA). OSA is primarily intended to be a federating platform for the simulation community: it is designed to favour the integration of new or existing contributions at every level of its architecture. The platform core supports discrete-event simulation engine(s) built on top of the ObjectWeb Consortium's Fractal component model. In OSA, the systems to be simulated are modeled and instrumented using Fractal components. Fractal components offer many advanced and original features, such as multi-programming language support and the ability to share sub-components. In OSA, the event handling is mostly hidden in the controller part of the components, which alleviates noticeably the modeling process, but also ease the replacement of any part of the simulation engine. Apart the simulation engine, OSA aims at integrating useful tools for modeling, developing, experimenting and analysing simulations. For this purpose it relies on the Eclipse development platform and its ability to be extended.

I. Introduction

For modeling and simulation of discrete event systems, the component approach was introduced in the 70's, with the DEVS formalism [START_REF] Zeigler | Theory of Modelling and Simulation[END_REF]. Indeed, an interesting property of the DEVS approach is that it allows to model complex systems by dividing the initial system, recursively and hierarchically, into smaller sub-systems. Since the DEVS formalism is a powerful way of formally describing complex systems, it has been widely used in discrete-event simulation softwares. However, DEVS do not address all the critical issues raised by component approaches. For example, as noted in preface of [START_REF] Zeigler | Theory of Modeling and Simulation[END_REF], the following critical issues for simulation have not received enough attention yet: model credibility (validation, verification, . . . ) and inter-operation (repositories, reuse of components, and resolution matching).

It is worth stressing that considering these major issues to be specific to the simulation domain would be too restrictive, because similar or identical issues are also being addressed in general purpose software [START_REF] Lau | A Taxonomy of Software Component Models[END_REF]. Indeed, the ability of components to ease the process of sharing and reusing code parts is clearly an attractive property for simulation models. However, this is also a common software issue that is addressed by most general purpose component models like CORBA [START_REF] Omg | [END_REF], COM [START_REF] Sessions | COM and DCOM: Microsoft Vision for Distributed Objects[END_REF] and Java-Beans [START_REF] Englander | Developing Java Beans[END_REF].

Consequently, building a new discrete-event simulation software based on a general purpose component architecture could noticeably help in solving the pending critical issues in the simulation domain. However, selecting the most appropriate general purpose component model to build such a new simulation software becomes a new issue. In this paper, we present our ongoing work on such a new simulation software, the Open Simulation Architecture (OSA). OSA relies on one of the latest emerging component models, the ObjectWeb's Fractal component model [START_REF] Bruneton | Recursive and dynamic software composition with sharing[END_REF], [START_REF] Bruneton | The fractal component model specification[END_REF].

The selection of the Fractal component model is no more than a heuristic answer to the previous selection issue. Indeed, the main design guideline of OSA is to pick as much as possible of the recent advances, emerging technologies, and trends from the general software engineering solutions and to apply them and study their benefits in the discrete-event simulation context. According to this policy, the Fractal model exhibits several appealing innovative features:

• it implements the separation of concerns paradigm [START_REF] Aksit | Separation and composition of concerns in the object-oriented model[END_REF];

• it allows the sharing of a sub-component between several distinct components (the same instance of a component may be inserted inside several distinct surrounding components);

• it offers multi-programming language support: for example, components written in the Java language may be mixed (or replaced) with components written in the C++ or C languages;

• if offers advanced mechanisms to support dynamic models (creation, destruction, and mutation of components), like factories components, template components, and binding primitives;

• it offers an extensible Architecture Description Language (ADL) based on the XML language.

The component model is a key part of OSA, but it is not enough to build a complete architecture. In addition to the Fractal component model, OSA adopts the following other elements and design principles:

• adopt both the Eclipse environment and Eclipse philosophy [START_REF] Des Riviêres | Eclipse: A platform for integrating development tools[END_REF]: provide an extensible environment open to contributions;

• extend the separation of concern principle to the whole architecture, identifying distinct user activities and strongly enforcing a separation between these activities;

• provide a model packaging and versionning system, and a repository service;

• Support distributed execution and middle-ware interconnexion (especially HLA [START_REF]IEEE Standard for Modeling and Simulation (M&S) High Level Architecture (HLA), Federate Interface Specification[END_REF]).

Before going further in the details of OSA, let us first emphasize that most of the previous elements are not yet implemented or even clearly specified. Indeed, because OSA is open to contributions, the OSA developments shall not (and will not) be conducted by a unique development team. In order to initiate the collaborative work process, the author's development team currently focus on a few basic elements, such as a simple simulation engine and tools for components edition and assembly.

In section II we first summarize the key features of the Fractal component model used in the OSA architecture. Then, through a simple modeling case study, section III illustrates how the Fractal model helps to enforce the previous objectives. Eventually, we describe the overall architecture of OSA in section IV and some of its key internals in section V.

II. The Fractal Component Model

From a software engineering point of view, a component is usually defined as "an independent unit of software deployment that satisfies a set of behavior rules and implements standard component interfaces that allow it to be composed with other components." [START_REF] Armstrong | Toward a common component architecture for high-performance scientific computing[END_REF] Fractal is the ObjectWeb Consortium component reference model [START_REF] Bruneton | Recursive and dynamic software composition with sharing[END_REF], [START_REF] Bruneton | The fractal component model specification[END_REF]. Fractal is neither a software environment nor a runtime executive. It is a specification. In other words, it is a set of rules and features that a component-based software architecture is supposed to follow or implement in order to be compliant with this model. Fractal does not mandate the use of any specific programming language. On the contrary, it allows to combine component implementations possibly based on different programming languages.

The Fractal specification defines several levels and sublevels of compliance. These levels allow an implementation not willing or not able to implement completely the model to state how much of the specification it complies with. At the lowest level, a component architecture claiming to be compliant with level 0.0 is just supposed to implement its components using the object programming paradigm. At the highest level, a component architecture claiming to be compliant with level 3.3 is supposed to fully implement all the features of the specification.

Hereafter, we summarize some of these key features (see [START_REF] Bruneton | The fractal component model specification[END_REF] for the complete description).

Component external structure. A Fractal component is an object-oriented unit of code that has external interfaces. These interfaces may be of two kinds: either client or server. The former emits service requests, the latter receives service requests. Interfaces are named. Their name must be unique for a given component but names may be reused for naming interfaces in other components. A client interface is intended to be bound to a server interface.

Hierarchical structure. Components may have a hierarchical structure (fig 1). Hierarchical components are made of a controller part (also called membrane) and a content part. The content part is composed of one or more components. Since a membrane and its content recursively form a component it may have external interfaces. It may also have internal interfaces. As external interfaces, internal interfaces may be either of type client, or of type server. Internal interfaces are only available to components of the content part. A component of the inner part may only bind its external interfaces to external interfaces of other inner components or to the inner interface of its surrounding controller. Therefore, the model strictly forbids a component to bind its external interfaces to the ones of components outside its membrane or inside its neighbouring (inner part) components. Interface Introspection. Introspection is the ability for an object to collect useful information about other objects (possibly including itself). In the Fractal model, components have the ability to introspect their interfaces. For example, a component may retrieve its own list of available internal and external interfaces.

Functional and controller interfaces. A functional interface is an interface used to offer or obtain services to or from other components. A controller interface is a server-only interface. It is offered to a component to access non-functional services, such as introspection, (re)configuration, persistence, service policy, life cycle control (ability to start/stop a component), and so on.

Factories and templates. A factory component is a component that has the ability to create other components. Fractal distinguishes two kinds of factories: generic factories, that have the ability to create several kinds of components, and standard component factories, that only have the ability to create one kind of component. Templates components are a special kind of standard factory components that may be recursively composed of factories, and serve as a model to create normal components in a quasi isomorphic manner (isomorphic meaning the created component has the same hierarchical structure as its creator template). Since factories are components and components are created from factories, a special component is required to initiate the recursion. This special component is a generic component factory called "bootstrap". Shared components. The Fractal model allows a component to appear in the content of several distinct enclosing components. Such components are called shared components. This property has two noticeable consequences: (i) a component is possibly placed under the control of several surrounding controller components and (ii) a shared component may directly interact with components located in the inner parts of several distinct components.

III. A simple modeling case study

This case study serves two objectives: (i) provide an example to illustrate the OSA modeling concepts, and (ii) demonstrate the usefulness of shared components.

In section III-A we first give a conceptual model of the system under study. Then, in section III-B, we present and discuss possible component implementations of this model with and without shared components.

A. Conceptual model

Figure 2 gives a conceptual view of our example system using a very simple form of Data Flow Diagram (DFD) [START_REF] Yourdon | Modern Structured Analysis[END_REF] in which we just show data flow interactions between system entities. The simple system is composed of two identical nodes that communicate with each other using the message passing paradigm. Each node of the system is decomposed in two sub-systems: the application and the transport. The application sub-system, itself, is decomposed in two sub-systems: the source and the sink.

The transport sub-system implements a reliable, ordered, connectionless, datagram routing and delivery service. It supports variable but limited datagram sizes.

The application sub-system contains a source subsystem that produces and sends new packets, and a sink sub-system that consumes and discards the packets received.

The source sub-system simulates file transfers. It repetitively executes the following actions: 1. sleep for a random time 2. pick a random file length flen; 3. Call ⌈ f len MAX P KT SIZE ⌉ times the transport's sending service to send fragments of at most MAX PKT SIZE bytes (the maximum packet size supported by the transport).

B. Component implementations

In the following, we discuss two particular implementation policies of our system example: a first one that does not use shared components and a second one that does use them.

B.1 Implementation without shared components.

The first implementation that comes to mind without shared components is the one shown on figure 3. This diagram reflects the client/server interactions between components which explains why some arrows have a reversed orientation compared to the ones of the data-flow diagram. For example, in the conceptual model the sink component receives data while in the client/server interaction diagram it initiates the reception service call.

The main quality of this implementation is to be structurally very close to the conceptual model. This good property is achieved by applying a very simple strategy in order to reflect interactions that need to cross subsystems boundaries: replicate the client and server interfaces of the interacting components on the internal and external sides of any surrounding membrane that need to be crossed. In this example, the consequence of the "replicate interfaces" strategy is that a node component have to expose an external transport interface. Therefore reusing node components in another system model would require a minimal understanding of this transport interface which contradict the separation of concerns principle: nodes cannot be used as self-contained "black-boxes" that hide their internal implementation details.

Furthermore, it is worth stressing that the simple "replicate-interface" strategy has several other negative effects that contradict the very fundamental philosophy of component-based design. For example, let's consider the evolve-ability good property of the component approach. The component-based approach is expected to ease the replacement of a component (when new improved versions of the component are released, for instance). For example, let us consider the following possible evolution of the node component: a new subcomponent modeling the physical network layer is added as shown in figure 4. Since interactions between transport sub-systems have been replaced by interactions between physical layer components, the component implementation of the new conceptual model implies modifications of all the surrounding components.

B.2 Implementation with shared components.

A possible implementation with a shared compoent is depicted in figure 5. In this new construct, the transport component is built hierarchically and contains two sub-components: the one named proxy is shared and the other, named local part is not. Indeed, the transport shall As a conclusion about shared components, philosophically, one may wonder whether using shared components is wise since it is a way of bypassing and somehow violating the components boundaries. It is worth stressing that in the context of simulation modeling, the isolation property of components may be considered from two point of views, which may be somehow confusing. The first point of view is the software engineering one and the second the simulation one. In [START_REF] Bruneton | The fractal component model specification[END_REF], the authors claim that "paradoxically shared components are useful to preserve component encapsulation". This is a software engineering point of view and exactly what we demonstrated with the simple model case study.

Let us now consider the simulation point of view. Sup-porting the shared component modelling feature is clearly breaking the tree structure of components. As stated in the Fractal specification, the component structure with shared components may be a directed acyclic graph and bindings (interaction paths) may form cycles. Since interactions are eventually translated into events, shared components may imply nasty causal effects, that should be properly handled by the discrete event simulation engine. However, while supporting this feature may add some complexity, especially for the engine, it also opens interesting perspectives, both for the engine implementation and for modelling.

The previous modelling case study gave an example of a modelling pattern in which the use of shared components may prove to be usefull. In section V, several places in the simulation engine implementation where the use of shared components may also prove to be usefull are discussed. Prior to entering into the details of the simulation engine implementation, an overall description of the OSA architecture is first given in the following section.

IV. The OSA architecture Studying a system using discrete-event computer simulations imply several activities [START_REF] Banks | Handbook of Simulation Principles, Methodology, Advances, Applications, and Practice[END_REF], [START_REF] Zeigler | Theory of Modeling and Simulation[END_REF], [START_REF] Jain | The Art of Computer Systems Performance Analysis: Techniques for Experimental Design, Measurement, Simulation, and Modeling[END_REF]. OSA aims at supporting a large number of these activities, which currently include (but are not restricted to) the one depicted on figure 6. A typical simulation study life-cycle is made of several iterations of the three following steps: 1. the simulation preparation, which involves tasks related to nearly all the functional concerns depicted on figure 6 except the analysis one; 2. the simulation execution (which belongs to the deployment functional concern); 3. the simulation run(s) analysis.

Notice that during the early iterations, while the model of the system is still under development, these three steps exist in a slightly different form: the simulation preparation consists in developing new models and building test scenarios, while the simulation run(s) analysis mainly consists in checking the validity of the models by comparing the effective outputs with the expected ones.

At any iteration of this life-cycle, the OSA architecture aims at providing a strong support for the two first steps, by means of an integrated (graphical) user interface and a simulation engine. Support for the last step, analysis, may be envisaged but is not a high priority because offthe-shelf power-tools already exist for this purpose.

The OSA software architecture follows a N-tier design which may involve the following layers: a front-end user interface layer, a middle-ware executive layer, a simulation engine layer, a data layer for results storage, and a source code management layer for models storage and versionning.

The three first layers are described further in the following sections.

A. Front-end user interface

The OSA architecture must provide tools to assist users in many tasks. Furthermore, the architecture should enforce a strong cooperation of these tools, using an integrated and easily extensible environment. For this purpose, we selected the Eclipse platform [START_REF] Holzner | Eclipse[END_REF], [START_REF] Des Riviêres | Eclipse: A platform for integrating development tools[END_REF].

Eclipse already provides a large amount of plugs-in to assist developers in various software development tasks: specification, development in several programming languages, unit testing, debugging, source code management, and so on. Some of these plugs-in are dedicated to the development of new Eclipse plugs-in, which explains the ever growing list of available plugs-in, and consequently its ever growing popularity.

An interesting feature of the Eclipse plugs-in is their ability to be extended. Indeed the Eclipse plug-in API defines extensions points [START_REF] Gamma | Contributing to Eclipse: principles, patterns, and plugs-in. The Eclipse series[END_REF]: plugs-in that implement such extension points (this is not mandatory) may be extended in order to build new enriched or specialized versions of the initial plugs-in.

Eclipse plugs-in are mainly used to build new Eclipse perspectives. An Eclipse perspective is dedicated organization of the Eclipse Graphical User Interface (GUI), offering support and specialized tools for a particular task.

Therefore, the development of the OSA user interface mainly consist in providing new Eclipse plugs-in and perspectives to support users in (possibly) all the tasks of the simulation study life-cycle. Within Eclipse, new plugs-in and perspectives may be developed concurrently, which enforces our collaborative development philosophy. Nevertheless, once new plugs-in and perspectives are made available for the community, this is up to each user to decide whether to install them or not.

B. Middle-ware layer

A middle-ware layer may optionally be used to support the execution of the simulations. Indeed, such architectures are often criticized for their potentially poor performance. Since performance is a critical issue for simulation, this architectural choice may be unwise. Therefore, the middle-ware layer is not mandatory in the OSA architecture, thanks to the Fractal component model: the distribution of component executions across a network through a middle-ware is an optional feature of Fractal that may, or may not be activated, without any change in the component functional implementations (the part of components that implements models).

The OSA architecture may support the distribution of the simulation executions across a network in different manners:

• distribution of several simulation-runs, each one executing on a single computer node. In this case, the distribution support required is very limited (a "gangscheduler" facility);

• distribution of one (or several) simulation runs across the network, simply using the Fractal model ability to distribute transparently the execution of the components, but without any cooperation of the simulation engine. Since the minimal requirement of the simulation engine, whether it executes in parallel or not, is to ensure consistency of event processing between components, this implements de facto a conservative mode of parallel execution [START_REF] Fujimoto | Parallel and distributed simulation systems[END_REF];

• distribution of one (or several) simulation runs across the network, using the Fractal model ability to distribute the execution of the components, and the cooperation of the simulation engine. Provided the components have the persistency non-functional feature in order to regularly save their global simulation state, this may lead to the optimistic modes of parallel execution [START_REF] Fujimoto | Parallel and distributed simulation systems[END_REF];

• the last form of distribution, which is somehow complementary of the previous ones, is achieved when the middle-ware is used to bridge together several simulation architectures, using the HLA standard for example [START_REF]IEEE Standard for Modeling and Simulation (M&S) High Level Architecture (HLA), Federate Interface Specification[END_REF].

C. Back-end simulation engine

The simulation engine is distributed over all the components that have a surrounding membrane implementing the simulation non-functional services. These services are accessed through a dedicated simulation-controller interface. Component with a simulation-controller interface are called managed components and those without such an interface are called passive components. This interface provides a unified access to a simulation-controller implementation. The simulation-controller implementation is not fixed. It may be totally or partially replaced.

The OSA simulation engine supports the three following semantics of interaction between components:

• synchronous interactions: the service requested by the client is synchronous with the simulation time (no simulated time consumption);

• asynchronous interactions: the client thread is not blocked while the service is being processed, concurrently, by the server. This asynchronous mode of interaction implies the client shall not expect a meaningful return value from its service call;

• blocking interactions: the client thread is blocked until the service is completed (or aborted) by the server. The client may expect a meaningful return value from its service call.

When at least one of the involved peers is a passive component, the only interaction mode available is the synchronous one.

The OSA simulation engine automatically translates asynchronous and blocking interactions into discreteevents, provided that both peers are managed components and thus each have a surrounding membrane with a simulation-controller interface.

The technique used to translate service invocations into simulation events was first introduced in [START_REF] Mussi | The prosit sequential simulator: a test-bed for object oriented discrete event simulation[END_REF]. It consists in using an object-oriented construction called a functor. A functor is the transformation of a method invocation into an object. Indeed, when a component issues a (functional) service request, through its client interface, its call is automatically intercepted by the controller part of the source component and reified into an event object.

This event object contains the current time in simulation, the event type, an object encapsulating the method called on the server side and its arguments (functor), and possibly other data that do not need to be further described.

In order to ensure a minimal inter-operability between components, the following rules are imposed :

• The particular event created at the time the service call is issued by the client is initialized with a type value of SOC (Start Of Call) and inserted in the event queue of the server interface;

• In case of a blocking service the event created at the time the service is completed by the server is initialized with a type value of EOR (End Of Reply) and inserted in the event queue of the client interface.

V. Implementing a simulation engine in OSA

As already mentioned in previous section, the simulation engine implementation is located (distributed) in the simulation-controller implementations that lay in each managed component.

However, because Fractal supports shared components, and controllers may be implemented as Fractal components, the simulation-controller, and thus the whole simulation engine may be shared between all managed components. In this case, the simulation engine becomes fully centralized.

Implementing a new simulation engine for OSA mainly consist in developing new event management policies and replacing the corresponding parts of the simulation-controller implementation.

The simulation-controller implementation is build using three types of abstractions, each having their own interface: the event queues, the event schedulers, and the explicit simulation services.

Dynamic structure of models [START_REF] Barros | Modeling Formalism for Dynamic Structure Systems[END_REF] may be implemented using the special factory and template components of the Fractal component model. Mutation of models and multi-models [START_REF] Ören | Dynamic Templates and Semantic Rules for Simulation Advisors and Certifiers[END_REF] may be implemented by replacing the implementation of the non-functional binding primitives of the components. This replacement is easy because the binding primitives are grouped together in a binding controller that lays in the controller-part of each component.

A. Event queues

An event queue should be associated with each functional interface of a component. But since event queues are components and components may be shared, all the interfaces of a component may share the same event queue. Or some may share a common queue and others may have their own.

B. Scheduling policy

Each managed component of the simulation should be associated with a scheduler. Like queues, schedulers may be shared amongst several components. The scheduler manage the current execution state of the component. A typical scheduler would support the following states :

• INIT; the component is in its initial state and awaiting for its start() service to be called;

• IDLE: the component is sleeping and waiting for an incoming request;

• RUNNING: the component is consuming processing time because it is servicing a request;

• BLOCKED: the component has issued a blocking service call and is waiting for the EOR event (service completion).

A typical blocking interaction between a client and a server component is shown on the temporal diagram of figure 7. As soon as several schedulers act concurrently, a synchronizing policy has to be implemented in order to enforce simulation time consistency [START_REF] Fujimoto | Parallel and distributed simulation systems[END_REF].

C. Explicit simulation services

During the simulation, the functional services (the service that implement the conceptual model of the system) may use the following services offered by the nonfunctional simulation-controller interface of their surrounding component:

• current time(): returns the current simulated time;

• terminate()/abort(): requests normal/abnormal termination of the simulation execution;

• wait until(time):

• wait until(condition):

• wait until(time,condition): the current component thread requests to stay in RUNNING state until the simulation time reaches time or the condition condition is matched. New incoming service calls may still be accepted by the component if the component has multiprocessing capability and an idle processing thread is available. Otherwise incoming calls may be either queued or refused according to the parameters of the queue of server interface on which the call was received.

• release(condition): resume execution of a component thread waiting for condition condition to be matched. This call is synchronous (non blocking): it returns true in case a thread is woken up and false otherwise.

• release(condition,time): resume execution of a component thread waiting for condition condition to be matched. This call is blocking its caller thread at most until time is reached. It returns true in case a thread is woken up and false otherwise (timeout time is reached). spin lock()/spin unlock(): primitive intended for components with multi-processing capability to ensure mutual exclusion between processing threads of the component (spinlock() blocks its caller processing thread as long as any other processing thread in the same component is not in the IDLE state).

VI. Conclusion

In this paper we introduced the Open Simulation Architecture and explained how this novel architecture could meet the expectations of a large part of the discrete-event simulation community. Indeed, the OSA architecture aims at (i) providing an open platform that supports simulationists in a wide set of their simulation activities, (ii) allowing the reuse and sharing of system models by means of a flexible component model, (iii) favouring contributions from the simulation community at any level of the architecture.

The critical issue of integrating numerous software contributions into a single open platform is addressed first by relying on a powerful and extensible GUI, Eclipse, whose philosophy is clearly to ease and favour such contributions and secondly by choosing the Fractal component model, which offers many interesting features, such as shared components and the ability to support multiprogramming language. Furthermore, the Fractal model adopts the separation of concerns paradigm, which ease the process of replacing any part of the architecture, from the model components to the very core parts of the simulation engine. Therefore, OSA may be used both as a testbed to experiment new parallel and distributed simulation engine algorithms or as an integration platform to reuse existing discrete-event simulation models.

It shall be noted that the OSA simulation engine does not enforce the use of any particular discrete-event formalism because the event management is kept hidden in the components non functional part. The resulting model of interactions is a classical (client-server) procedural programming model, which may serve as a basis to support or encapsulate most popular modeling formalisms, such as Petri nets or DEVS.

In order to initiate third party contributions in OSA, a detailed specification of the simulation-controller interfaces and a first implementation of the simulation engine, based on the Java language and Aspect Oriented Programming techniques, are planned to be re-leased during the first quarter of 2006. This implementation is based on the AOKell implementation of the Fractal model [START_REF] Seinturier | Recent developments in AOKell[END_REF].

VII. Acknowledgements

This work is co-sponsored by the french National Research Agency and the INRIA Research Institute.

Fig. 1 :

 1 Fig. 1: Example of Fractal hierarchical component.

Fig. 2 :

 2 Fig. 2: Simple system conceptual model.

Fig. 3 :

 3 Fig. 3: Component implementation without shared component

Fig. 4 :

 4 Fig. 4: A possible evolution of the simple system

Fig. 5 :

 5 Fig. 5: Implementation with a shared component

Fig. 6 :

 6 Fig. 6: OSA functional architecture.

Fig. 7 :

 7 Fig. 7: Temporal diagram of interactions during a blocking service call.


Lihan Xie
email: lihan.xie@insa-lyon.fr

S Ébastien Baguet
email: sebastien.baguet@insa-lyon.fr

Benoit Prabel
email: benoit.prabel@cea.fr

R Égis Dufour
email: regis.dufour@insa-lyon.fr

Parametric analysis of the nonlinear behavior of rotating structures

In nonlinear rotordynamics, techniques can take advantage of the periodic steady state behavior to predict quickly and accurately the mass unbalance response to a series of parameters, especially with the presence of certain nonlinearities which leads to nonlinear dynamics and complicated responses. The method proposed here calculates the response curve by combining Harmonic Balance Method (HBM), Alternating Frequency-Time (AFT) method and continuation. The singular points where a stability change often arises are detected with the sign change of the Jacobian determinant and then located through a penalty method that increases the solving equation system by a completing constraint. Tracking these points, which provides an efficient way to analyze parametrically the nonlinear behavior of a system, can be fulfilled, once again, by the continuation technique.

INTRODUCTION

Generally speaking, the rotating systems utilized in the energy production have a small rotor-stator gap, are able to run during long periods, and are mounted on hydrodynamic * Address all correspondence for other issues to this author. bearings. Rotor-stator interactions in case of blade loss, crack propagation due to fatigue, and a variable stiffness due to the nonlinear restoring forces of the bearings can make the rotordynamics nonlinear and the responses complicated: significant amplitude and frequency shifts are introduced, suband super-harmonics appear, and hysteresis occurs. It is of great importance to understand, predict and control this complicated dynamics. This is commonly achieved by means of numerical simulations due to its efficiency, low cost, repeatability, operability, etc. Reaseachers have investigated numerically the dynamic of systems with various nonlinear properties, such as rotor-stator contact [START_REF] Groll | The harmonic balance method with arc-length continuation in rotor/stator contact problems[END_REF] [2] [START_REF] Peletan | Quasi-periodic harmonic balance method for rubbing self-induced vibrations in rotor tator dynamics[END_REF], crack breathing [START_REF] Chen | Nonlinear response and dynamic stability of a cracked rotor[END_REF] [START_REF] Cheng | The influence of crack breathing and imbalance orientation angle on the characteristics of the critical speed of a cracked rotor[END_REF], hydrodynamic bearings [START_REF] Khonsari | Stability boundary of non-linear orbits within clearance circle of journal bearings[END_REF] [START_REF] Brown | Chaos in the unbalance response of journal bearings[END_REF], on board rotor mounted on hydrodynamic bearings [START_REF] Dakel | Nonlinear dynamics of a support-excited flexible rotor with hydrodynamic journal bearings[END_REF], etc.

At the design stage of nonlinear mechanical systems, a particular attention must be paid to the influence of parameters in order to reach the more efficiently as possible the configuration which leads to optimal runs. Moreover, appropriate parameters choices can avoid dangerous resonance phenomena by moving the resonance frequencies out of the operating frequency range or by decreasing the resonance amplitudes to acceptable levels. Therefore, direct paramet-2 EQUILIBRIUM PATH ric analysis which quickly assesses the effect of parameters on the dynamic behavior is necessary for numerical investigation of nonlinear systems.

In vibration analysis, first of all, the resonance levels and frequencies are of primary interests in the majority of investigations of forced response. Many works aim at finding the peak vibration amplitudes in nonlinear dynamical systems. Petrov [START_REF] Petrov | Direct parametric analysis of resonance regimes for nonlinear vibrations of bladed disks[END_REF] [START_REF] Petrov | Method for sensitivity analysis of resonance forced response of bladed disks with nonlinear contact interfaces[END_REF] applies the multiharmonic balance method to compute the worst vibration cases of bladed disks with friction contact interfaces. In [START_REF] Liao | Maximization of the vibration amplitude and bifurcation analysis of nonlinear systems using the constrained optimization shooting method[END_REF], Liao combines the shooting method and Floquet theory along with a Global Search algorithm to determine the resonant peak of nonlinear systems. Similarly, Balaram et al [START_REF] Balaram | Optimal design of multi-parametric nonlinear systems using a parametric continuation based genetic algorithm approach[END_REF] have developed a normminimizing strategy in the coupling of continuation and Genetic Algorithm to perform an optimal design.

Then, determining the local stability of a periodic solution is particularly interesting in an engineering context since only stable solutions are experimentally encountered [START_REF] Nayfeh | Nonlinear oscillations[END_REF]. Moreover, a change in the stability can lead to significant, qualitative, and possibly dramatic changes in the system response. Efforts have been previously made by researchers for accurate detection of the singular points which include the limit (turning) points and bifurcation points. Since the tangent stiffness matrix becomes singular at these points, most of the proposed methods use this characteristic as the critical constraint which is added to the governing equations to form a so-called augmented system defining singular points [START_REF] Moore | The calculation of turning points of nonlinear equations[END_REF] [15] [START_REF] Battini | Improved minimal augmentation procedure for the direct computation of critical points[END_REF]. Filippa proposes regularization penalty of increased system to avoid the singularity of the augmented system when approaching a bifurcation point [START_REF] Felippa | Traversing critical points by penalty springs[END_REF].

This work presents a rapid systematic method for the direct parametric analysis of the targeted vibration and the change of stability of nonlinear systems. First, the nonlinear dynamic analysis is performed in the frequency domain owing to the high computational efficiency of harmonic balance method (HBM) [START_REF] Groll | The harmonic balance method with arc-length continuation in rotor/stator contact problems[END_REF]. With the help of continuation method, all dynamic equilibrium solutions of nonlinear systems are determined. Then, for stability analysis, the Floquet theory is employed [START_REF] Liaw | Dynamic stability and chaos of system with piecewise linear stiffness[END_REF] [19] [START_REF] Peletan | A comparison of stability computational methods for periodic solution of nonlinear problems with application to rotordynamics[END_REF]. The determinant of the Jacobian matrix is monitored because it helps to pick out the bifurcation points from limit points [START_REF] Felippa | Traversing critical points by penalty springs[END_REF]. The singular points are then located by adding a new constraint equation to the solving system.

Next, the equation of forced response levels and bifurcation points can be calculated directly as a function of nonlinear parameters or excitation level thanks to the added constraint function and, once again, the continuation method.

EQUILIBRIUM PATH

The equation of motion for the forced vibrations of a structure with nonlinear characteristic takes the form:

M ẍ(t) + Cẋ(t) + Kx(t) + f nl (x, ẋ) = p(ω,t) (1)
where x(t) is a vector of displacements for all n degrees of freedom (DOF), M, C, K are generalized n × n mass, damping, and stiffness matrices of finite element model, p(ω,t) is a vector of excitation forces, which, in our case, is periodic (unbalance force for rotating machines), and ω is the excitation frequency. The rotation effects, such as gyroscopic moments, Coriolis forces, centrifuge effect, etc, can be taken into account in the damping and stiffness term. Harmonic balance method (HBM) is utilized for its efficiency with respect to time domain method. The time variation of displacements and external forces for the steady-state periodic regime are represented by restricted Fourier series:

x(t) = X 0 + N ∑ k=0 X k c cos(kωt) + X k s sin(kωt) (2)
p(t) = P 0 + N ∑ k=0 P k c cos(kωt) + P k s sin(kωt) (3)
where

X = [X 0 , X c 1 , X s 1 , ..X c N , X s N ] T , P = [P 0 , P 1 c , P 1 s , ..P N c , P N s ]
T are the vectors of Fourier coefficients for displacements and external excitation with truncated decomposition up to order N.

The main idea of harmonic balance method is to identify the periodic solution through estimations of the truncated Fourier coefficients by solving the nonlinear equation of motion in frequency domain. By substituting (2) (3) into (1), then applying a Galerkin procedure, the nonlinear differential equation ( 1) is transformed into a nonlinear algebraic equation system of dimension n × (2N + 1):

R(X, ω, λ) = Z(ω)X + F NL (X, λ F ) -P(ω, λ p ) = 0 (4) where Z = diag(K, Z 1 , ..Z k , ..Z N ) with Z k = K -k 2 ω 2 M ωC -ωC K -k 2 ω 2 M
, F NL and P are vectors of harmonic coefficients for nonlinear forces and excitation forces. λ F and λ p are parameters of nonlinear and excitation forces that can be varied later for parametric analysis. A robust scheme, which obtains a set of discrete points to trace the equilibrium path in the range of parameters characterizing the structure, is based on Newton's method. The utilized Newton-Raphson method consists in an incrementaliterative procedure for approaching the exact solution

JδX k = -R(X k ) (5)
X k+1 = X k + δX k (6)
where J(X, ω, λ) = ∂R k ∂X = Z(ω) + ∂F NL ∂X is the Jacobian matrix, R(X k ) the residue of the equilibrium equation of motion at iteration k.

PARAMETRIC CONTINUATION OF SINGULAR POINTS -STABILITY BOUNDARY DETERMINATION

The Fourier coefficients of nonlinear efforts F NL and its derivative ∂F NL ∂X involved in the Newton-Raphson iterations are obtained by Alternating Frequency-Time (AFT) method [START_REF] Cameron | An alternating frequency/time domain method for calculating the steady-state response of nonlinear dynamic systems[END_REF]. The AFT scheme uses fast direct and inverse Fourier transforms to compute the nonlinear forces in the time domain and then switch back to the frequency domain for the reason that nonlinear forces are usually much easier to evaluate in the time domain than in the frequency domain.

Nonlinear systems often have several possible responses for a given excitation frequency ω.

The pseudo-arc length continuation method [START_REF] Crisfield | A fast incremental/iterative solution procedure that handles snap-through[END_REF] combined with the abovementioned algorithm permits to follow the solution branch beyond turning points in order to obtain both stable and unstable solutions of response curves. The prediction step of continuation is firstly performed in the tangent direction to the solution curve. Then, this estimation is corrected iteratively in the orthogonal direction for the purpose of cancelling the residue R 1 :

R 1 (X, ω) = J c δX k+1 δω k+1 = J ∂R k ∂ω ∆X T i ∆ω i δX k+1 δω k+1 = -R k 0 (7)
where [∆X T i ∆ω i ] and δX k+1 δω k+1 t stands for the unitary tangent vector and the orthogonal correction vector respectively, k and k + 1 correspond to the iteration number. The pseudo-arc length continuation method thus adds an equation and an unknown to the solving system. The iterative calculation is accomplished when the required accuracy is achieved. Consequently, the equilibrium path on the study range is obtained.

LOCATION OF SINGULAR POINTS AND STA-BILITY ASSESSMENT

Singular points that often accompanied by a change of stability, are generally distinguished into turning (limit) points and bifurcation points. In some cases, turning points are located near the vibration resonance. Besides, the analysis of bifurcation points has large practical importance since any of the multitude of branching solutions can be realized for the same operating conditions.

The determinant values of Jacobians J and J c provide the information to classify the found solutions as follow:

regular point i f det(J) = 0 and ∂R T ∂ω y = 0 turning point i f det(J) = 0 and ∂R T ∂ω y = 0 bi f urcation point i f det(J) = 0 and ∂R T ∂ω y = 0 ( 8
)
where y is the eigenvector associated with the zero eigenvalue of J. During the continuation of the response curve, a sign change of det(J c ) between two consecutive points indicates the presence of a bifurcation point.

To locate singular points precisely, it is more efficient to characterize the singular points using the equation Jy = 0 rather than the determinant of J. Between the two consecutive points, the one which has a smaller absolute value of det(J) is chosen to approach the exact singular point. which leads to solve the following augmented system for each iteration of Newton-Raphson:

  J 0 ∂R ∂ω ∂(Jy) ∂X J ∂(Jy) ∂ω 0 2y T 0      δX δy δω    = -    R + γe j Jy y 2 -1    (9) with J = Z + ∂F NL ∂X ∂(Jy) ∂X = ∂ 2 F NL ∂X 2 y ∂R ∂ω = ∂Z ∂ω X - ∂P ∂w ∂(Jy) ∂ω = ∂Z ∂ω y
For locating a bifurcation points, an additional equation

∂R T
∂ω y = 0 and an new unknown γ are added. The term γe j in ( 9) is added only when locating a bifurcation point where e j is a unit vector with j-th component equal to 1, γ is an auxiliary variable for eliminating the singularity of Jacobian [START_REF] Felippa | Traversing critical points by penalty springs[END_REF] and its initial value is assumed to be 0. The approximation of the eigenvector y is calculated by performing singular value decomposition of the Jacobian and selecting a vector corresponding to the minimum magnitude of the singular value.

This algorithm applies only to simple singular points. In more complex cases, such as the change of dynamic regime form periodic to quasi-periodic (see §5.2), the detection of the bifurcation point can be performed by means of a classical stability assessment.

Assessment of the local stability consists in applying a small perturbation to the equilibrium solution, and then checking whether or not this perturbation subsides with time. Floquet theory is the most widely used method for evaluating the stability of solution. There are several algorithms, both in the time and frequency domain, to compute the eigenvalues of the so-called monodromy matrix [START_REF] Peletan | A comparison of stability computational methods for periodic solution of nonlinear problems with application to rotordynamics[END_REF].

In this work, the monodromy matrix Φ = [φ φ] is obtained by integrating the following system:

M φ(t) + (C + ∂f nl ∂ ẋ ) φ(t) + (K + ∂f nl ∂x )φ(t) = ∂p ∂x (10)
over one period of motion with the Newmark algorithm.

The eigenvalues of the monodromy matrix, also known as the Floquet multipliers, provide information on the stability of the periodic solution. If all the multipliers are inside the unit circle in the complex plane (R e , I m ), the corresponding solution is asymptotically stable and called a stable limit cycle. If at least one multiplier is outside the unit circle, there is a local bifurcation of limit cycle with a loss of stability.

PARAMETRIC CONTINUATION OF SINGULAR

POINTS -STABILITY BOUNDARY DETERMINA-TION Once the singular point (X 0 , ω 0 ) is specifically detected along the equilibrium branch, it is possible to follow its vari-5.1 Duffing oscillator

NUMERICAL RESULTS

ation with respect to a parameter λ(λ F , λ p ) with a continuation method similar to that described in paragraph §3. In order to do this, λ is considered as a new unknown and an equation r 2 (X, ω, λ) = 0 is added to the augmented system [START_REF] Petrov | Direct parametric analysis of resonance regimes for nonlinear vibrations of bladed disks[END_REF], which requires the Newton-Raphson corrections to be orthogonal to the prediction step. The continuation of the singular point is thus conducted by successive corrections in the following form:

       J 0 ∂R ∂ω e j ∂R ∂λ ∂(Jy) ∂X J ∂(Jy) ∂ω 0 ∂(Jy) ∂λ 0 2y T 0 0 0 0 ∂R T ∂ω 0 0 0 ∆X T 0 ∆ω 0 ∆λ                   δX δy δω δγ δλ            = -            R + γe j Jy y 2 -1 ∂R T ∂ω y r 2            (11) with ∂ R ∂λ = ∂F NL ∂λ - ∂P ∂λ ∂(Jy) ∂λ = ∂ 2 F NL ∂X∂λ y (12)
Again, the term γe j and δγ are necessary to eliminate the singularity of the augmented system during the bifurcation point continuation.

NUMERICAL RESULTS

Numerical simulations are presented in this section in order to demonstrate the validity and the advantage of the proposed approach.

Duffing oscillator

The Duffing oscillator is a one-degree of freedom system with a nonlinear restoring force proportional to the cube of the displacement. The accuracy of the proposed method in predicting the periodic motions of the forced Duffing oscillator is examined by comparing the stable branch obtained by HBM (19 harmonics) with those obtained from direct temporal integration. The system is described by the differential equation:

m ẍ(t) + c ẋ + kx(t) + αx(t) 3 = p 0 cos(ωt) ( 13
)
where α is the nonlinear stiffness coefficient, and p 0 is the excitation amplitude.

Using the numerical continuation of periodic solutions with the excitation frequency as the continuation parameter, the assessment is firstly carried out with the parameters: m = 1, c = 0.1, k = 1, α = 1, p 0 = 0.5. The frequency response curve is shown in Fig. 1. A classical resonance curve can be observed while the blue dots and red dots stand for stable and unstable solutions respectively. The two magenta points are obtained by monitoring the sign change of the conditions in [START_REF] Dakel | Nonlinear dynamics of a support-excited flexible rotor with hydrodynamic journal bearings[END_REF]. Since only det(J) passed through 0, these are limit points with stability change and no bifurcation occurs. By applying the continuation technique presented above, the two limit points are followed when the nonlinear coefficient α is varied from 0.01 to 10. In Fig. 2, the response curve of the linear system (α = 0) is plotted in black line, the blue dotted curve nearby is the system response for α 0 = 0.01 (slightly nonlinear) which is totally stable, while the response on the right is calculated for α 1 = 10 (strongly nonlinear). The green dotted-line represents the direct tracking of the first limit point while the magenta dotted-line tracks the lower limit point. The zone that between these two dotted-lines is the unstable area of the Duffing system when its nonlinear coefficient varies. The two lines are the stability boundaries defining amplitude jump: when frequency increases, am amplitude jump arises when the system vibration level comes close to the green line, and when frequency decreases, an amplitude jump will be predicted by the red line.

When the nonlinear effect is stronger, or the excitation force is higher, bifurcation points are found as shown by green circles in Fig. 3. Besides of the two limit points previously studied, there are other limit points at the superharmonic resonances.

Nonlinear Jeffcott rotor

The second test case is a modified Jeffcott rotor which can come into contact with a stator that is modeled as an added stiffness [START_REF] Groll | The harmonic balance method with arc-length continuation in rotor/stator contact problems[END_REF] [START_REF] Peletan | Quasi-periodic harmonic balance method for rubbing self-induced vibrations in rotor tator dynamics[END_REF]. The rotor is made of a weightless shaft carrying a disk with mass m at the middle of the span. The clearance between the rotor and the stator is denoted by h. The stator, which is rigidly fixed, has an elastic contact surface modeled as a symmetrical set of radial springs with isotropic stiffness, k c . The equations of motion are shown below

m ẍ + c ẋ + kx + k c (1 - h r )(x -µysign(v rel )) = p b ω 2 cosωt m ÿ + c ẏ + ky + k c (1 - h r )(µxsign(v rel ) + y) = p b ω 2 cosωt ( 14
)
where k is the stiffness of the shaft, r = x 2 + y 2 is the radial displacement, p b is the imbalance amplitude and v rel is the relative velocity between the rotor and the stator at the contact point.

All the calculations are carried out with the same conditions: m = 1, c = 5, k = 100, k c = 2500, h = 0.105, p b = 0.1, except that the friction coefficient µ is considered as the var- ied parameter. Firstly, the response curve of the system when µ = 0 is plotted in Fig. 4. Both the detection of sign change of Jacobian determinant and the stability assessment were performed. As shown, the red dots indicate the unstable branch while the magenta points mark the limit points. The response has the same shape for µ = 0.1.

Yet, when the friction coefficient increases to 0.12 then to 0.2, the stability distribution changes a lot. In Fig. 5, the singular points are marked, however, the stability change next to the dimensionless frequency 0.3 doesn't cause a sign change in the Jacobian determinant. In fact, it's related to the regime change of the system from synchronous full annular rub motion to partial rub motion, which means the stable periodic motion changes to quasi-periodic motion [START_REF] Jiang | Determination of the global responses characteristics ofapiecewise smooth dynamical system with contact[END_REF]. This fact indicates that monitoring the sign change of Jacobian is efficient to provide information on simple singular points only (saddle-node, asymmetric, and pitchfork bifurcations), but it is insufficient to predict all the stability change of system solutions. Constraint equations and associated extended systems for tracking complex bifurcations (period doubling, Neimark-Sacker) can be found in [START_REF] Doedel | Computation of periodic solution bifurcations in odes using bordered systems[END_REF]. Their use in the framework of the HBM are beyond the scope of this paper. Figure .6 shows another stable branch between the two unsta- [START_REF] Peletan | A comparison of stability computational methods for periodic solution of nonlinear problems with application to rotordynamics[END_REF] ble branches indicating that partial rub motion changes back to synchronous full annular rub.

Figure .7 collects all the calculated responses for several values of friction coefficient µ. Blue and red lines indicate stable and unstable branches respectively. The half-green half-magenta line is calculated with the continuation technique which permets a direct following of the bifurcation points with respect to µ. The stability analysis of converged solution on this line distinguishes stable solutions (green) from unstable ones (magenta). As observed, the continuation technique succeeds to locate directly all the bifurcation points. Where they correspond to a change of stability in the response curve, they are stable (green branch), whereas they are unstable (magenta curve) when there is no stability change involved.

Finite element nonlinear rotor

A finite element nonlinear rotor model [START_REF] Lalanne | Rotordynamics prediction in engineering[END_REF] [20] was also tested (see Fig. 8). This rotor is composed of 13 Timoshenko beam elements. Two linear isotropic mountings are located at both ends of the rotor (k yy = k zz = 6e7N.m -3 , c yy = c zz = 600N.s.m), and three disks (modeled by additional mass and inertia) are added to nodes 3, 6 and 11. The aim is to calculate its response with an imbalance force located on disk 2 (m b .g = 2e -2kg.m). The rotor makes contact with a circular, static, and rigid stator located in the vicinity of disk 2. The contact model is identical to that de-

Mẍ(t) + C(ω)ẋ(t) + Kx(t) + f c (x) = p(ω,t) (15)
where M, C(ω) and K represent the mass, damping and stiffness matrices respectively. C(ω) includes the gyroscopic matrix, which varies with ω.

The calculations for finite element rotor are carried out in Cast3m [START_REF]CAST3[END_REF]. Figure.9 represents the response curve with stable and unstable branches for the friction coefficient µ = 0. The black line stands for linear response. Two turning points are also detected. The response curve for µ = 0.03 is shown in Fig. 10. Same as in Fig. 5, the system has passed stable and unstable phases due to complex motion changes due to the strong non-linearity arising from contact. Similarly, not all bifurcations have been detected.

In Fig. 11 the magenta line is directly calculated with the continuation technique and covers all the bifurcation points with respect to µ (from 0 to 0.2). This example illustrates that the continuation technique is able to locate directly the bifurcation points in applications of finite element analysis.

CONCLUSIONS

An efficient frequency-domain method is presented for the rapid parametric analysis of stability changes in nonlinear dynamic systems, which provides directly the stability boundary with respect to parameters such as the system nonlinearity or excitation level.

First, the equilibrium solutions are calculated by continuation technique combined with harmonic balance method. Then stability of equilibrium solutions is determined by Floquet theory. The determinant of Jacobian provides information on limit points and bifurcation points. After locating these singular points, the proposed method is able to follow them directly with respect to a varied parameter. The application examples have illustrated the efficiency and accuracy of proposed method.

Numerical developments are fulfilled in both Matlab and Cast3m, paving the way for application of the method to the nonlinear dynamics of rotors modelized with 3D finite elements.

Fig. 1 .Fig. 2 .Fig. 3 .

 123 Fig. 1. Frequency-response curve of the Duffing system for m = 1, c = 0.1, k = 1, α = 1, p 0 = 0.5. (blue dot : stable solution; red dot : unstable solution; magenta point : turning point)

Fig. 4 .

 4 Fig. 4. Forced response curve of Jeffcott rotor for µ = 0. (blue dot : stable solution; red dot : unstable solution; magenta point : turning point; green circle: saddle-node bifurcation)

Fig. 5 .Fig. 6 .

 56 Fig. 5. Forced response curve of Jeffcott rotor for µ = 0.2. (blue dot : stable solution; red dot : unstable solution; magenta point : turning point; green circle: saddle-node bifurcation)

5. 3 Fig. 7 .Fig. 8 .

 378 Fig. 7. Singular point tracking of Jeffcott rotor as a function of friction coefficient µ. (blue branch : stable solution; red branch : unstable solution; green point : stable solution; magenta point : unstable solution )

Fig. 9 .

 9 Fig. 9. Forced response curve of finit element rotor for µ = 0 (blue dot : stable solution; red dot : unstable solution; magenta point : turning point)

Fig. 10 .

 10 Fig. 10. Forced response curve of finite element rotor for µ = 0.03 (blue dot : stable solution; red dot : unstable solution; magenta point : turning point)

Fig. 11 .

 11 Fig. 11. Singular point tracking of finite element rotor as a function of friction coefficient µ. (blue branch : stable solution; red branch : unstable solution; magenta point : singular points )