B Enaux

V Granet

O Vermorel
email: olivier.vermorel@cerfacs.fr

C Lacour

L Thobois

Large Eddy Simulation of a motored single-cylinder piston engine: numerical strategies and validation

Keywords: Internal combustion engine, Multi-cycle Large-Eddy Simulation, Cycle-to-cycle variations, Unstructured grids, Acoustic

This paper describes a compressible Large Eddy Simulation (LES) used to investigate cyclic variations for nonreacting flow in an optical single cylinder engine setup. The simulated operating point is part of a large experimental database designed to validate LES for cycle-to-cycle prediction, and constitutes a first step towards the realization of fired operating points. The computational domain covers almost the whole experimental setup (intake and exhaust plenums, intake and exhaust ducts, cylinder) to account for acoustic phenomena. The assessment of the computation is performed in two regions of the domain: the intake and exhaust duct predictions are compared to the results of a Helmholtz solver and the experiment (pressure transducers and Particle Image Velocimetry (PIV)) while the in-cylinder dynamics are confronted to PIV measurements. The ability of the developed methodology to capture the correct level of cycle-to-cycle variations

Introduction

A major challenge for the development of Internal Combustion (IC) engines is to improve fuel economy and to reduce pollutant emissions while maintaining or enhancing engine performances. New strategies using stratified charge direct injection or downsizing with high levels of exhaust gas recirculation have this potential, but can impact on the combustion stability and trigger high Cycle-to-Cycle Variations (CCV). Nowadays Computational Fluid Dynamic (CFD) applied to engine development mainly relies on the Reynolds Averaged Navier Stokes (RANS) approach which solves ensemble or phase-averaged mean quantities. However, this method does not provide information about individual cycles. The study of cycle-to-cycle variability is thus impossible and other numerical methods such as LES (Large Eddy Simulation) must be employed. LES of laboratory-scale but also real turbulent combustors, such as gas turbine chambers, have recently appeared in the literature [START_REF] Moin | Large-eddy simulation of realistic gas turbine combustors[END_REF][START_REF] Boileau | LES of an ignition sequence in a gas turbine engine[END_REF][START_REF] Patel | Simulation of spray-turbulence-flame interactions in a lean direct injection combustor[END_REF][START_REF] Fureby | Towards the use of large eddy simulation in engineering[END_REF]. In comparison to numerous publications related to gas turbines for example, LES of IC engines is still a field to explore. Two major reviews, one by Haworth [START_REF] Haworth | Large-eddy simulation of in-cylinder flows[END_REF] and the other by Celik et al. [START_REF] Celik | Large eddy simulations of in-cylinder turbulence for internal combustion engines: a review[END_REF], summarize the work accomplished by the LES engine community. They both highlight the potential of LES for predicting the details of the transient in-cylinder flow dynamics and demonstrate that LES is able to produce superior predictions than any of the current RANS approaches. They also support the view that LES can provide essential information in regard to the nature and origins of CCV. However, they both underline that computational cost is a major shortcoming that LES will have to overcome to ensure its practical use. Indeed, by nature, CCV can only be studied by performing multiple cycles and statistical analysis. [START_REF] Lax | Difference schemes for hyperbolic equations with high order of accuracy[END_REF] cycles are reported to be a minimum for a correct sampling in terms of mean values and 50 cycles are probably necessary to yield reliable statistics in terms of fluctuations [START_REF] Heywood | Internal combustion engine fundamentals. McGraw and Hill Series in Mechanical Engineering[END_REF][START_REF] Celik | Assessment measures for URANS/DES/LES: an overview with applications[END_REF]. Fortunately, the constant increase in computational resources and the use of massively parallel machines now make such computations realistic. Several multi-cycle LES-based computations have been reported in the literature since the review papers of Haworth [START_REF] Haworth | Large-eddy simulation of in-cylinder flows[END_REF] and Celik et al. [START_REF] Celik | Large eddy simulations of in-cylinder turbulence for internal combustion engines: a review[END_REF]. Hasse et al. [START_REF] Hasse | Numerical investigation of cyclic variations in gasoline engines using a hybrid URANS/LES modeling approach[END_REF] performed 13 consecutive cycles of a motored engine using Detached Eddy Simulation (DES). Cycle-to-cycle velocity fluctuations at the spark plug were investigated at the time of ignition. It was shown that the differences in velocity could lead to significant combustion variations when computing the corresponding reacting case. Considering 5 different initial scenarios, Goryntsev et al. [START_REF] Goryntsev | Large eddy simulation of fuel-air mixing in a direct injection SI engine[END_REF][START_REF] Goryntsev | Large eddy simulation based analysis of the effects of cycle-to-cycle variations on air-fuel mixing in realistic DISI IC-engines[END_REF] performed 10 consecutive nonreacting two-phase flow LES cycles for each scenario. It was found that CCV phenomena were directly linked to turbulence and that velocity cyclic variations had a great impact on the fuel-air mixing process. Multicycle simulations of fired operating points have also been performed recently [START_REF] Richard | Towards large eddy simulation of combustion in spark ignition engines[END_REF][START_REF] Thobois | Using LES to investigate reacting flow physics in engine design process[END_REF][START_REF] Vermorel | Towards the understanding of cyclic variability in a spark ignited engine using multi-cycle LES[END_REF]. In particular, it was shown in [START_REF] Vermorel | Towards the understanding of cyclic variability in a spark ignited engine using multi-cycle LES[END_REF] that LES was able to reproduce the experimental trend in terms of CCV, even with only 9 simulated cycles. In this study, the sources of CCV were investigated and the importance of the turbulence generated by the tumble breakdown was highlighted. However, most authors highlight the lack of experimental data specifically dedicated to the validation of LES [START_REF] Celik | Large eddy simulations of in-cylinder turbulence for internal combustion engines: a review[END_REF][START_REF] Dugue | Applicability of large eddy simulation to the fluid mechanics in a real engine configuration by means of an industrial code[END_REF][START_REF] Vermorel | Towards the understanding of cyclic variability in a spark ignited engine using multi-cycle LES[END_REF]. The French project SGEmac has been specially built to validate LES of CCV by acquiring a large experimental database on different operating points of a four-valve spark ignition engine [START_REF] Lacour | Exploring cyclic variability in a spark-ignition engine using experimental techniques, system simulation and large-eddy simulation[END_REF]. This database includes motored engine and fired operating points, with low or high levels of CCV. A detailed characterization of the various operating points is achieved with numerous temperature and temporal pressure measurements at several locations of the set-up. Optical diagnostics (Particle Image Velocimetry (PIV), OH Laser Induced Fluorescence (LIF), combustion chemiluminescence) allow also the investigation of the local properties of the flow. The computational domain includes most of the experimental setup: cylinder, intake and exhaust ducts, intake and exhaust plenums. This strategy has potentially two major drawbacks, its CPU cost since the computational domain is large (∼ 2.5 m long) and its modeling cost since some additional devices (flame-arrestors for example) must be accounted for. But it also offers two great advantages: first, it allows to capture acoustic phenomena in the ducts which are a possible source of CCV [START_REF] Haworth | Large-eddy simulation of in-cylinder flows[END_REF][START_REF] Dugue | Applicability of large eddy simulation to the fluid mechanics in a real engine configuration by means of an industrial code[END_REF][START_REF] Vermorel | Towards the understanding of cyclic variability in a spark ignited engine using multi-cycle LES[END_REF]; second, it enables a proper definition of intake and exhaust boundary conditions. Indeed, in most studies, boundary conditions are imposed in the intake and exhaust ducts, close to the cylinder [START_REF] Goryntsev | Large eddy simulation of fuel-air mixing in a direct injection SI engine[END_REF][START_REF] Richard | Towards large eddy simulation of combustion in spark ignition engines[END_REF][START_REF] Vermorel | Towards the understanding of cyclic variability in a spark ignited engine using multi-cycle LES[END_REF]. With such a compact domain, several challenging issues need to be tackled such as turbulence injection, back-flow treatment or acoustic wave reflection for instance. With an extended domain, these issues vanish naturally. In this paper, only a motored engine case is reported. It constitutes an essential validation step before computing more complex operating points with combustion. This operating point is experimentally known to be very stable with low CCV. Therefore, the objective here is not only to discuss cyclic variations but also to prove that the chosen methodology is suitable to describe with high fidelity the engine behavior, from an aerodynamic and acoustic point of view. To ensure a proper statistical sampling of the flow field, 27 consecutive LES cycles are performed. This paper is organized as follows. The experimental setup and engine settings are recalled in section 2. The LES solver and the models are presented in section 3. Section 4 describes the meshing methodology, the prescription of the boundary conditions and the modeling of the flame-arrestors. Section 5 presents a comparison of the LES results with PIV acquisitions and temporal pressure measurements specifically introduced in the experimental setup. The flow behavior in the intake and exhaust lines is first analyzed before considering the in-cylinder predictions.

Experimental setup

The SGEmac engine test bench operated at IFP has been specifically designed to investigate cyclic combustion variability. This single-cylinder is a four-valve, pent-roof, spark ignition engine with a flat piston. The whole setup is displayed in Fig. 1(a) and the main parameters are summarized in Table 1. The operating points acquired on this bench are documented in [START_REF] Lacour | Exploring cyclic variability in a spark-ignition engine using experimental techniques, system simulation and large-eddy simulation[END_REF]. For fired points, the test bench is fueled with gaseous propane while for the non-reactive point the engine is only filled with air. Air and propane flow rates are controlled by sonic nozzles. Air is introduced in a tranquilization plenum and propane is then added in a mixing plenum. At the engine exhaust, gases are tranquilized in a third plenum. Since the engine is operated in a fully premixed mode, a flame-arrestor is added for safety reasons before each plenum (Fig. 1(a)). One hundred experimental consecutive cycles were acquired for each operating point. Figure 1(b) shows the location of the various experimental acquisitions. PIV measurements are performed through a quartz cylinder and an optically transparent intake duct (planes 5 and 6) providing instantaneous and phase-averaged velocity fields. Several pressure transducers are used to record at one CAD resolution the pressure evolution along the intake (noted 1, 2 and 3) and exhaust ducts (noted 4).

In the literature, different parameters are used to quantify CCV. Variations in Indicated Mean Effective Pressure (IMEP) or in-cylinder peak pressure P max are the most usual indicators of CCV [START_REF] Ozdor | Cyclic variability in spark ignition engines[END_REF][START_REF] Young | Cyclic dispersion in the homogeneous charge spark-ignition engine. a literature survey[END_REF]. For the present motored engine case (no combustion), IMEP is not a relevant parameter and thus the variations can only be characterized with the in-cylinder flow field and P max . Experimental measurements have shown that this operating point is very stable in terms of P max with relative variations of about 1%.

LES of piston engines flows

LES is performed with a fully compressible Navier-Stokes solver on unstructured grids with a cell-vertex finite-volume formulation [START_REF] Moureau | Numerical methods for unsteady compressible multi-component reacting flows on fixed and moving grids[END_REF][START_REF] Schmitt | Large-eddy simulation and experimental study of heat transfer, nitric oxide emissions and combustion instability in a swirled turbulent high-pressure burner[END_REF][START_REF] Gourdain | High performance parallel computing of flows in complex geometries -part 1: methods[END_REF][START_REF] Gourdain | High performance parallel computing of flows in complex geometries -part 2: applications[END_REF]. Centered spatial schemes and explicit time-advancement are used to control numerical dissipation and capture acoustics. For the present case, a Lax-Wendroff scheme [START_REF] Lax | Difference schemes for hyperbolic equations with high order of accuracy[END_REF], 2nd order in space and time, is employed with a time step conditioned by the acoustic CFL number. The LES solver takes into account changes of heat capacity with temperature and composition using tabulated values of heat capacity for each species. The multispecies fluid follows the ideal gas law, P = ρrT where P is the pressure, ρ is the density, T is the temperature and r is the mixture gas constant. Viscous transport terms are handled using classical gradient approaches. The fluid viscosity and the heat diffusion coefficient follow Sutherland's law and Fourier's law respectively. The species diffusion coefficients are obtained using a species Schmidt number along with Hirschfelder and Curtis approximation and velocity corrections for mass conservation [START_REF] Poinsot | Theoretical and Numerical Combustion[END_REF]. Filtering the Navier-Stokes equations yields unresolved terms, which need to be modeled. The following closures of the subgrid scale (SGS) terms are used in the LES solver:

-The SGS stress tensor τ i, j t is closed using a Boussinesq approximation. This approximation introduces a turbulent viscosity ν T which is estimated with the standard Smagorinsky model [START_REF] Smagorinsky | General circulation experiments with the primitive equations: 1. the basic experiment[END_REF] with a constant C S = 0.18.

-The SGS species fluxes J k i t are modeled using a species SGS turbulent diffusivity D t k = ν T /Sc t k , where Sc t k is the turbulent Schmidt number (Sc t k =0.6 for all species k). -The SGS energy flux q i t is modeled using the SGS turbulent thermal conductivity λ t th = ρν T C p /Pr t with C p the heat capacity at constant pressure of the mixture and Pr t the turbulent Prandtl number (Pr t = 0.6).

The methodology of Cook and Cabot [START_REF] Cook | Hyperviscosity for shock-turbulence interactions[END_REF] is used to treat shocks by introducing a hyper-viscosity in the viscous stress tensor which thickens the shock front. Such shocks can appear at moments where exhaust valves are opened. The boundary condition treatment is based on a multi-species extension [START_REF] Moureau | Numerical methods for unsteady compressible multi-component reacting flows on fixed and moving grids[END_REF] of the Navier-Stokes Characteristic Boundary Conditions (NSCBC) [START_REF] Poinsot | Boundary conditions for direct simulations of compressible viscous flows[END_REF].

Numerical Setup

Acoustics are a potential source of CCV and must be captured by LES. In the project SGEmac, two methodologies are explored [START_REF] Lacour | Exploring cyclic variability in a spark-ignition engine using experimental techniques, system simulation and large-eddy simulation[END_REF]. The first considers the coupling of a system simulation (1D) and LES (3D), where the 3D simulation is only applied on the combustion chamber and a part of the intake and exhaust pipes. The present scope of this paper concerns the second methodology which consists in applying LES on most of the experimental setup. Figure 1(a) shows the computational domain. The computed intake line extends from the mixing plenum to the cylinder while the exhaust line is fully computed. The two flame-arrestors are accounted for. The air plenum and the line between the two intake plenums have been removed since it is expected that most of the acoustic activity takes place between the cylinder and the mixing plenum. This configuration allows a simple and unambiguous definition of the boundary conditions (two plenum conditions). The following subsections detail the key points of the methodology used in this work.

Mesh management

The whole computation domain is meshed using tetrahedra: it covers about 2.5 m from the intake boundary condition to the exhaust boundary condition and includes two plenums of 20 l (intake) and 50 l (exhaust). In the present study, an engine cycle is divided into 41 phases (i.e. 41 grids) in order to maintain the quality of the mesh throughout the cycle. Between two computation phases, a second-order interpolation is used. The grids of each phase are generated before the simulation using a two-staged smoothing approach (Laplacian and optimization-based smoothing) [START_REF] Forsythe | Validation of a fluid-structure interaction model for a bileaflet mechanical heart valve[END_REF]. As a reference, recent and comparable LES of IC engine used 20 phases with hexahedral-based grids [START_REF] Richard | Towards large eddy simulation of combustion in spark ignition engines[END_REF][START_REF] Vermorel | Towards the understanding of cyclic variability in a spark ignited engine using multi-cycle LES[END_REF] and 35 phases with hybrid tetrahedral/prismatic-based grids [START_REF] Hasse | Numerical investigation of cyclic variations in gasoline engines using a hybrid URANS/LES modeling approach[END_REF]. This method allows to avoid distortions of tetrahedra and to preserve accuracy and stability [START_REF] Vermorel | Towards the understanding of cyclic variability in a spark ignited engine using multi-cycle LES[END_REF]. The smallest mesh contains 2.2 million cells at top-dead-center (TDC) while 9.6 million cells are required for the biggest mesh at bottom-dead-center (BDC). The characteristic cell size into the cylinder is of the order of 0.8 mm. The meshes are specifically refined around the valve seats, with a minimum of seven cells within a residual valve lift of approximatively 0.35 mm. From the intake/exhaust ports to the plenums, the meshes are gradually derefined with a biggest characteristic cell size in the plenums of about 20 mm. Note that meshes are not refined at the walls since a law-of-the-wall formulation is used [START_REF] Schmitt | Large-eddy simulation and experimental study of heat transfer, nitric oxide emissions and combustion instability in a swirled turbulent high-pressure burner[END_REF]. Figure 2 shows an overview of a typical tetrahedral mesh during the intake stroke.

Boundary conditions

At the top of the intake and exhaust plenums, a pressure imposed boundary condition is set. The target pressures are not constant but time-varying plenum signals extracted from the experimental data which change between 4.35 × 10 4 and 4.5 × 10 4 Pa for the inlet plenum and 9.9 × 10 4 to 1.019 × 10 5 Pa for the outlet plenum. These signals are quasi sinusoidal with a frequency dictated by the engine rotational speed. As mentioned in the introduction, usual IC engine simulations consider boundary conditions directly in the ducts, close to the cylinder. Compared to a plenum signal, signals measured in the ducts (i.e. closer to the cylinder) vary much more and contain a wider spectrum of scales with many high frequencies: including the exhaust and intake ducts and imposing pressure in the plenums is obviously making the simulation easier. Now, even if a larger computational domain makes the simulation simpler, it does not solve all boundary condition problems because acoustic impedances at inlets and outlets must be controlled: a section where pressure is strictly imposed is fully reflective for acoustic waves propagating in the LES while it is not necessarily so in the experiment. There is no easy solution to this difficulty. The best one is to locate boundary conditions as far as possible from the engine in places where pressure is really constant. Here, quasi nonreflective NSCBC are used via a Linear Relaxation Method (LRM) [START_REF] Rudy | A non-reflecting outflow boundary condition for subsonic Navier Stokes calculations[END_REF][START_REF] Poinsot | Boundary conditions for direct simulations of compressible viscous flows[END_REF] to impose pressure in the two plenums. LRM boundary conditions are well known to act as first order low-pass filters introducing a cut-off frequency that separates waves that will be reflected from the ones that will leave the domain [START_REF] Selle | The actual impedance of non-reflecting boundary conditions: implications for the computation of resonators[END_REF]. This cut-off frequency is directly proportional to the relaxation coefficient imposed on the target pressure signal. Thus, one of the main advantages of a plenum boundary condition compared to a duct boundary condition is to allow the use of lower relaxation coefficients since the imposed signal varies much less and contains lower frequencies. As a consequence, the plenum boundary condition will be acoustically less reflective. However, the introduction of a low relaxation coefficient can deteriorate the target signal, both in magnitude and in phase. In other words, the pressure signal which is actually recovered in the computation may differ from the target signal imposed in the plenums. For the range of relaxation coefficients used in this study, preliminary tests have shown that the LRM induces a negligible damping in magnitude. However, it leads to a phase shift between the target pressure (measured experimentally) and the inlet pressure obtained in the LES. This problem can be solved by shifting the target pressure in time as follows: for a signal composed of a single frequency f , Selle et al. [START_REF] Selle | The actual impedance of non-reflecting boundary conditions: implications for the computation of resonators[END_REF] have demonstrated that the phase shift Φ introduced by the reflexion coefficient A of a characteristic boundary condition with a relaxation coefficient K can be estimated by:

Φ = arctan Im(A) Re(A) with A = - 1 1 -i 2π f K (1)
Knowing Φ, it is then possible to deduce the corresponding temporal shift ∆t to impose on the initial target signal in order to obtain, after application of the LRM, the right desired signal:

∆t = Φ/ (2π f ) (2)
A similar methodology was extended to the multi-frequency pressure signals of the intake and exhaust boundary conditions. Figure 3 presents the results for the intake boundary condition. In this case, the shift induced by the relaxation coefficient (Eqs. 1 and 2) is evaluated to -0.0417 s i.e. -50 CAD. The experimental pressure signal has then been shifted by -50 CAD and specified as target signal in the computation. The resulting inlet pressure in the LES, after applying boundary condition corrections, matches the experimental signal correctly. Thus this methodology allows to obtain acoustically well-defined boundary conditions while preventing any drift of the imposed signal. In addition, this boundary condition can take into account the acoustics in the ducts which is a possible source of CCV.

Modeling of the flame-arrestors

The flame-arrestors were introduced in the experiment for obvious security reasons, but from the computational point of view they represent an additional difficulty. Preliminary experimental studies of the flame-arrestor have been performed on an acoustic test bench [START_REF] Tran | Analysis and control of combustion instabilities by adaptive reflection coefficients[END_REF]. They have shown that acoustic in the ducts is strongly perturbed by the presence of these flame-arrestors indicating that they

∆ P = ρl ′ σ sin(α) ∂U W ∂ t + 1 2C D 2 ρ U W 2 σ 2 sin 2 (α) (3)
where l ′ = πd 4 +h is a characteristic length determined to match the Howe model [START_REF] Eldredge | Numerical investigation of the acoustic behavior of a multi-perforated liner[END_REF], U W is the bulk velocity and α is the inclination angle of the aperture. In this expression, the discharge coefficient C D depends on the bulk velocity U W (given by the manufacturer).

Results and Discussion

The whole simulation covers 27 consecutive cycles. Computation starts at IVO with a flow at rest. The turnaround time of a single LES cycle is about 18 hours on 400 processors of a SGI Altix ICE 8200 cluster. Figure 5 presents the evolution of the trapped mass for the 27 consecutive LES cycles. The first two cycles are strongly influenced by the initialization of the simulation and are discarded in the following analysis. Thus, 25 cycles remain for the comparison with the experimental data. The 25 remaining LES cycles are very stable in terms of trapped mass around 200.6 mg. The difference compared to the value estimated in the experiment (199.9 mg) is approximately 0.5%.

The validation of the multi-cycle LES is organized as follows: -First, the flow field activity in the ducts is validated using PIV measurements (plane 6 in Fig. 1(b)) and temporal pressure acquisitions (probes 1 and 4 in Fig. 1(b)) in the intake and exhaust manifolds. Acoustic wave propagation in the intake and exhaust ducts is specifically investigated. -Second, the LES quality inside the cylinder is assessed and the flow behavior is compared to the experimental data in terms of in-cylinder pressure and aerodynamic field throughout the cycle. CCV are highlighted and analyzed.

Flow behavior in the intake and exhaust ducts

Figure 6 shows the pressure evolution during the whole cycle at the intake and exhaust manifolds (probes 1 and 4 respectively in for the intake duct and between EVC and EVO for the exhaust duct, see Table 1).

The fact that LES and experimental signals match very well during the whole cycle shows that the modeling described in Section 4.3 is able to reproduce the right impedance of the flame-arrestors. The signal extracted from a simulation without flame-arrestor (also reported in Fig. 6 (right) for the exhaust line) shows that the main effect of the flame-arrestor is to damp the oscillations, the phase of the signal being little affected. The same behavior may be observed for the intake line. Despite this satisfactory global behavior, the LES signal in the intake duct shows a superimposed high frequency signal around 430 Hz which is not present in the measurements.

To determine the nature of the 430 Hz frequency appearing in the LES pressure signal, two diagnostics were used: a spectral analysis of the LES results via a spectral map and an acoustic analysis of the configuration via a Helmholtz solver. The spectral map consists in FFT analysis applied on multiple pressure probes recorded on the centerline of the intake duct. The sample frequency of the pressure signals is 100 kHz and the spectral resolution of the FFTs is 0.4 Hz since 25 LES cycles are recorded corresponding to a total physical time of 2.5 s. The spectral map of the LES pressure signals reveals a very wide spectrum with many characteristic frequencies, going (for the most energetic ones) from 10 Hz, the frequency of the valve opening, to 630 Hz. The structure of the 430 Hz mode is shown in Fig. 7 (top). This mode is an acoustic eigenmode of the intake duct and this can be checked by using an acoustic Helmholtz solver [START_REF] Nicoud | Acoustic modes in combustors with complex impedances and multidimensional active flames[END_REF]. This code (called AVSP) solves the eigenvalue problem associated to the Helmholtz equation which solves the wave equation in a flow:

∇.(c 2 ∇ p) + ω 2 p = 0 ( 4
)
where c is the local speed of sound and p is the Fourier transform of the pressure fluctuations at the radial frequency ω = 2π f . This solver has been applied to the intake line with the valves closed since AVSP is not able to take into account moving boundaries. The results in Table 2 show that four main longitudinal acoustic modes coexist in the intake line. The modes 1, 3 and 4 are respectively the quarter wave, three-quarter wave and five-quarter wave modes of the full line. The mode 2 is the quarter wave mode between the valves and the intake manifold and it is damped in the LES. The analysis of the spatial structure of the five-quarter wave mode at 405 Hz (Fig. 7 bottom) shows that the correspondence with the 430 Hz mode found in the LES is excellent. The slight shift of frequency and of spatial structure is certainly due to the fact that the spectral map is built over the whole cycle, i.e. valves opened and closed, contrary to the acoustic configuration. This result confirms that this 430 Hz LES mode is a real mode of acoustic nature and not a hydrodynamic effect or a numerical artifact. Of course, the AVSP result cannot explain why this acoustic mode is damped in the experiment and not in the computation. Experimental measurements made without flame-arrestors have shown that they are not the source of damping. This point must still be explored but the agreement between LES and experiment in Fig. 6 remains excellent anyway.

Figure 8 compares the phase-averaged crank-resolved pressure obtained in the LES to the experimental measurements at probe 3. The corresponding phaseaveraged axial velocity is also shown. Experimental velocity values are extracted from PIV measurements on the centerline of the duct. Negative velocities indicate a fluid flowing from the plenum to the cylinder. LES and experiments are in good agreement during the whole cycle, both for the pressure and the velocity. The strong back-flow occurring at IVO (CA = -360 • ) is clearly visible. This phenomenon is not localized close to the cylinder only. During about 20 CAD the complete intake line flows back to the plenum. Thereafter, the flow enters the cylinder until a second minor back-flow occurs again just before IVC (CA = -170 • ). When the valves are closed at CA = -120 • , the flow continues to oscillate between the intake ports and the plenum. The variations from one cycle to another are roughly constant during the whole cycle, of the order of ±2 m/s. It is worth reminding that this behavior corresponds to a probe located on the centerline of the duct. A probe located at the same longitudinal position but shifted in the radial direction exhibits noticeably different evolution. In extreme cases, both positive and negative velocities may be found in a given section of the duct as shown in Fig. 9. This complex evolution of the velocity, rapidly varying from negative to positive values both in space and in time, confirms how challenging it would be to set a boundary condition in section 3 (Fig. 1(a)) because it is too close to the engine. At this section, the flow changes sign with time and is not homogeneous in space. Considering Fig. 9, the simple question of defining if the boundary condition should be an inlet or an outlet is already a difficult task. Moreover, it has to be noted that all hot gases flowing back to the intake boundary condition are "lost" for the computation since the boundary condition is generally built only to feed cold fresh gases into the domain. This is not very problematic for a motored engine computation where all gases have a quasi constant composition and temperature but it may become an important source of errors for reacting cases where gases leaving the engine strongly differ from fresh gases feeding it.

In-cylinder predictions

Flow characteristics

Figure 10 shows instantaneous velocity fields obtained by LES at four distinct times. These fields exhibit the typical structures encountered during the intake, compression and exhaust strokes. At IVO, a back-flow first occurs due to the lower pressure within the intake pipe. Then, after about 20 CAD, fresh gases start entering the cylinder. The annular valve jets interact with the cylinder head and the piston to form smaller structures (Fig. 10(a)). At higher valve lifts (Fig. 10(b)), the jet direction coupled to the piston movement induces a well-known tumbling motion in the cylinder. During compression, the tumble motion becomes more and more structured until TDC as illustrated in Fig. 10(c). This organized motion disappears during the expansion stroke. The exhaust stroke also begins with a back-flow from the exhaust pipe into the cylinder. Due to the very small valve lift, the flow is sonic around the valve heads (Fig. 10(d)). At higher exhaust valve lifts, the pressures of the exhaust pipes and the cylinder balance out and air is driven good agreement. The LES behavior around TDC indicates that the compression ratio used in the simulation differs from the experimental one, certainly because of blowby issues often encountered experimentally with optical engines.

Assessment of quality and reliability of the multi-cycle LES

Before going into more details, the quality and the reliability of the multi-cycle LES has to be assessed. The quality of a LES may be investigated using two categories of estimators: single-computation estimators or grid-model procedures which require several simulations (typically a standard simulation, then a simulation with a modified SGS model, and finally a coarse grid simulation) [START_REF] Klein | An attempt to assess the quality of large eddy simulations in the context of implicit filtering[END_REF][START_REF] Celik | Index of resolution quality for large eddy simulations[END_REF][START_REF] Klein | Assessment of les quality measures using the error landscape approach[END_REF].

In the present study, the efforts in terms of grid meshing and CPU time are too important for considering the second class of indicators. Therefore the quality of the simulation is evaluated by means of the ratio between SGS viscosity ν T and fluid viscosity ν. most of the cycle. The most critical periods in terms of flow resolution are encountered around compression TDC and EVO. Indeed, during the compression stroke, the grid resolution stays almost constant (from 0.8 mm at BDC to 0.65 mm at TDC) while the turbulent scales strongly decrease (typically by a factor of ten [START_REF] Heywood | Internal combustion engine fundamentals. McGraw and Hill Series in Mechanical Engineering[END_REF][START_REF] Celik | An assessment of turbulence scales relevant to internal combustion engines[END_REF]), leading to an insufficient resolution around TDC. At EVO, the grid is much more refined (∼ 0.07 mm in the valve seats) but the turbulent activity generated by the high speed jets in the valve seats (∼ 480 m/s) is too strong to avoid higher levels of turbulent viscosity. The reliability of the simulation is also controlled by the statistical sampling.

In order to compare the LES to the experiments, phase averaging is performed: phase-averaged quantities are obtained using 25 cycles in the LES and 100 cycles in the experiment. Previous studies [START_REF] Moureau | Towards large eddy simulation in internal-combustion engines: simulation of a compressed tumble flow[END_REF][START_REF] Dugue | Applicability of large eddy simulation to the fluid mechanics in a real engine configuration by means of an industrial code[END_REF][START_REF] Toledo | Large eddy simulation of the generation and breakdown of a tumbling flow[END_REF][START_REF] Hasse | Detached eddy simulation of cyclic large scale fluctuations in a simplified engine setup[END_REF] have already addressed the impact of the number of cycles used for phase averaging. For gauging the relevance of the LES sample size, a convergence analysis was performed first on the experimental data. It consists in evaluating the difference between a phase-averaged velocity over n cycles (n = 10, 25, 50, 75) and the reference phase-averaged velocity over 100 cycles. The difference ε CA n evaluated by Eq. 5 is computed along a vertical velocity profile sketched in Fig. 16(a):

ε CA n = 1 ∆ z cylinder head piston U 2 x +U 2 z CA n - U 2 x +U 2 z CA 100 U 2 x +U 2 z CA 100 dz , (5)
where . CA n is the ensemble averaging over n cyles at a given crank angle CA and ∆ z is the integration vertical height. Table 3 shows the error ε CA n at four instants of the intake and compression strokes. It suggests that n = 25 is sufficient to evaluate the mean in-cylinder dynamics since the experiment shows a difference of about 6% between a phase-averaging over 25 cycles and a phase-averaging over 100 cycles. A similar study was carried out for the root mean square (RMS) of the velocity. Table 3 indicates that a minimum of 50 cycles would be necessary to obtain an error of 10% on RMS velocities, thus confirming the hypothesis of Heywood [START_REF] Heywood | Internal combustion engine fundamentals. McGraw and Hill Series in Mechanical Engineering[END_REF] or Celik et al. [START_REF] Celik | Assessment measures for URANS/DES/LES: an overview with applications[END_REF] for example. It is worth noting that 25 cycles allow an acceptable estimate of the RMS values nonetheless (∼ 10% error except in the early times after valve opening). The same analysis can be then done for LES results but only on 25 cycles. Table 4 shows that computing 10 cycles only (as done by several authors [START_REF] Dugue | Applicability of large eddy simulation to the fluid mechanics in a real engine configuration by means of an industrial code[END_REF][START_REF] Vermorel | Towards the understanding of cyclic variability in a spark ignited engine using multi-cycle LES[END_REF][START_REF] Hasse | Numerical investigation of cyclic variations in gasoline engines using a hybrid URANS/LES modeling approach[END_REF]) allows a reasonable analysis of the mean flow. Conversely and as expected, this is insufficient for RMS quantities considering that the reference values obtained with 25 cycles are already a limited estimation of the real values.

It has to be mentioned that this analysis is of course restricted to an aerodynamic point of view since the configuration is a motored engine case. For fired operating points, the same kind of analysis could be carried out with combustion indicators like the flame surface or heat release fields. If the same conclusions are probably valid for stable and weakly unstable operating points, one can expect a noticeably different conclusion for highly unstable cases. (i.e. a much larger statistical sampling would be necessary).

In-cylinder flow field

Cycle-to-cycle analysis

In this section, the emphasis is put on the CCV analysis of pressure and velocity variations in the cylinder: CCV are usually studied in terms of maximum cylinder pressure but this quantity gives only a global overview of the engine stability. The study of velocity CCV allows a more precise and local characterization of the flow unsteadiness.

Exhaust Intake

Pressure variations

Figure 15 presents the return map of the maximum in-cylinder pressure. This consists in plotting the maximum pressure of cycle i + 1 (normalized by the maximum pressure of the mean cycle) versus the maximum pressure of the preceding cycle i.

The level of CCV in the cylinder pressure is very low for both LES and experiment (relative variation < 0.4%), meaning that LES correctly reproduces the stability of the experimental point. One can observe that neither LES nor experiment exhibit a regular pattern meaning that no obvious relation exists between a cycle and the following one. However, little less CCV are observed in the LES compared to the experiments. Contrary to the experiment where the pressure signal at the plenum slightly varies from cycle to cycle, the same signal is imposed in the LES for all cycles. Tests are currently underway to check if this simplification could modify the results in terms of cyclic variations. The difference in sampling size could be another reason for these discrepancies. Finally, the use of a relatively low-order numerical scheme combined with a Smagorinsky model known to be dissipative may also lead to an excessive dissipation that could reduce CCV.

Aerodynamics variations

Even if pressure CCV are very low, CCV can be significant for the in-cylinder flow field. To illustrate this point, Fig. 16 shows CCV in the axial velocity at different crank angles. The axial velocity component is plotted along the z-axis in the center of the cylinder. For the sake of clarity, the 100 experimental cycles are not superimposed and only a confidence envelope is plotted. The lower and upper part of the confidence envelope contains 95% of the experimental cycles and are estimated as Ux 100 -2 • σ Ux and Ux 100 + 2 • σ Ux respectively, with σ Ux the standard deviation of the axial velocity. Fig 16(b) shows that CCV levels are slightly higher in the upper region of the profiles during the first times of the intake stroke, due to valve jets. During compression, Fig. 16(d) and (e) show that the variations are homogeneous over the entire height of the chamber with lower levels than during the intake stroke. However, a normalization by the mean velocity reveals that the relative intensity of the variations is stronger during the compression stroke, as already shown by Goryntsev et al. [START_REF] Goryntsev | Large eddy simulation of fuel-air mixing in a direct injection SI engine[END_REF]. The comparison between the computed instantaneous cycles and the confidence envelope is satisfactory during the whole cycle, especially during the intake stroke, meaning that LES is able to capture the variations of the incylinder dynamic. To determine precisely where aerodynamic CCV occur during a cycle, a proper orthogonal decomposition (POD) [START_REF] Lumley | Stochastic Tools in Turbulence[END_REF] analysis of the in-cylinder velocity field is carried out using the snapshot method introduced by Sirovich [START_REF] Sirovich | Turbulence and the dynamics of coherent structures[END_REF]. To take into account the deformation during a cycle, a phase-dependent POD is performed every 5 CAD [START_REF] Fogleman | Application of the proper orthogonal decomposition to datasets of internal combustion engine flows[END_REF]. Figure 17 presents the energy distribution of the ten highest energy modes found in the POD as a function of the CA. The largest part of the energy (∼ 70%) is contained in the first mode which corresponds to the mean flow (similar results were found in [START_REF] Baby | A new methodology to analyse cycle-to-cycle aerodynamic variations[END_REF]). The other modes reveal the velocity fluctuations and their corresponding energy denotes their level. In other words, at a considered CA, POD allows an energetic classification of the CCV of the in-cylinder dynamic. The ratio between energies of the first and second mode suggests that CCV are very low during the intake stroke (∼ 3%) and become stronger during the com- pression with a maximum value near TDC (∼ 9.5%) due to the vortex breakdown and near BDC of the expansion stroke (∼ 7.5%). This result is thus consistent with the previous analysis for velocity profiles. Figure 18 presents iso-surfaces of the most and second most energetic POD modes found inside the cylinder. It shows that CCV at 280 CAD before TDC are mainly created in the mid-plane of the cylinder by the shear between the valve jets. Then the flapping of the jets becomes the main source of variation as illustrated in Fig. 18(b).

At 180 CAD before TDC, the location of the most energetic mode completely changes. Indeed, at this crank angle the tumbling motion becomes the most energetic mode instead of the intake jets. During the compression stroke, Fig. 18(d) shows that variations are provided by the wall/tumble interaction and some variations are also concentrated in the center of the cylinder, which may be related to a tumble core precession [START_REF] Maurel | Extended proper orthogonal decomposition: Application to jet/vortex interaction[END_REF][START_REF] Borée | Disruption of a compressed vortex[END_REF].

Conclusion

This paper presents a methodology designed to study CCV (Cycle-to-Cycle Variations) in IC engines using multi-cycle LES. In order to limit the influence of the boundary conditions and to ensure proper acoustic resolution, the computational domain takes into account almost the whole experimental bench including intake and exhaust ducts and plenums. The numerical strategy includes mesh movement of tetrahedral elements, modeling of the flame-arrestors and specific treatments of boundary conditions at inlet and exhaust. The methodology is tested on a nonreacting stable operating point with low CCV and compared to experimental results. 27 consecutive LES cycles are computed but only 25 cycles are considered for the analysis. Very reasonable turnaround time is achieved since only 18 hours (elapsed time) are necessary to compute a single cycle. Acoustic activity in intake and exhaust ducts is validated against both pressure transducer measurements and a Helmholtz solver. The flow dynamic in the cylinder is well reproduced by LES with levels of CCV comparable to the experiments. The highest absolute variations are found during the intake stroke but a normalization by the mean flow reveals that the maximum relative variations rather occur during the compression stroke. This report is confirmed by a proper orthogonal decomposition analysis that spatially locates the dynamic variations and denotes their corresponding energy. The high levels of dynamic CCV have little impact on the evolution of the cylinder pressure. The variations in maximum pressure are indeed very low both in the LES and in the experiments. The stable nature of the operating point is thus well reproduced by the simulation. The next step is to apply this methodology for fired operating points with higher CCV levels. The ability of LES to predict both stable and unstable points and their corresponding CCV will be analyzed in future publications.

Fig. 1

 1 Fig. 1 (a) Sketch of the experimental SGEmac engine test bench and locations of the LES boundaries. (b) Location of the experimental measurements: (1,2,3,4) pressure transducers; (5) PIV measurement in a vertical cross-section of the combustion chamber; (6) PIV measurement in a vertical cross-section of one of the intake ports.

Fig. 2

 2 Fig. 2 View of a typical tetrahedral mesh during the intake stroke.

Fig. 3

 3 Fig. 3 Pressure signals at the intake boundary condition. Experimental signal ( ), shifted target signal imposed in the computation ( ), signal recovered after imposition of the quasi non reflective NSCBC ( ).

Fig. 4

 4 Fig. 4 Flame-arrestor: (a) in the experiment; (b) in the computation.

Fig. 5

 5 Fig.[START_REF] Celik | Large eddy simulations of in-cylinder turbulence for internal combustion engines: a review[END_REF] Evolution of the trapped mass in the cylinder for the 27 LES cycles. The experimental trapped mass ( ) is estimated to 199.9 mg.

  Figure6shows the pressure evolution during the whole cycle at the intake and exhaust manifolds (probes 1 and 4 respectively in Fig1(b)) and the good agreement between phase-averaged LES and experimental signals. The two signals exhibit the same progressive damping when the valves are closed (between IVC and IVO for the intake duct and between EVC and EVO for the exhaust duct, see Table1). The fact that LES and experimental signals match very well during the whole cycle shows that the modeling described in Section 4.3 is able to reproduce the right impedance of the flame-arrestors. The signal extracted from a simulation without flame-arrestor (also reported in Fig.6(right) for the exhaust line) shows that the main effect of the flame-arrestor is to damp the oscillations, the phase of the signal being little affected. The same behavior may be observed for the intake line. Despite this satisfactory global behavior, the LES signal in the intake duct shows a superimposed high frequency signal around 430 Hz which is not present in the measurements.To determine the nature of the 430 Hz frequency appearing in the LES pressure signal, two diagnostics were used: a spectral analysis of the LES results via a spectral map and an acoustic analysis of the configuration via a Helmholtz solver. The spectral map consists in FFT analysis applied on multiple pressure probes recorded on the centerline of the intake duct. The sample frequency of the pressure signals is 100 kHz and the spectral resolution of the FFTs is 0.4 Hz since 25 LES cycles are recorded corresponding to a total physical time of 2.5 s. The spectral map of the LES pressure signals reveals a very wide spectrum with many characteristic frequencies, going (for the most energetic ones) from 10 Hz, the

Fig. 6

 6 Fig. 6 Phase-averaged mean pressure in the intake (a) and exhaust (b) manifolds, measurement noted 1 and 4 respectively in Fig. 1(b). LES: ; experiments: . LES results without modeling of the flame-arrestors ( ) are also shown for the exhaust manifold.

Fig. 7

 7 Fig. 7 Spatial structure of the five-quarter wave longitudinal mode in the intake ducts. Comparison of the field of |p'| between spectral map on the duct axis (top) and the Helmholtz solver results in the whole intake domain (bottom).

Fig. 8

 8 Fig. 8 Phase-averaged crank-resolved pressure (a) and velocity (b) in the intake duct, measurements noted 3 and 6 respectively in Fig. 1(b). LES: ; experiments: .

Fig. 9

 9 Fig.9Phase-averaged axial velocity field (LES) in a cutting plane through the intake ducts at 180 CAD before TDC (location 3 in Fig.1(b)). The isoline of null velocity is also plotted.

Fig. 10

 10 Fig. 10 Instantaneous velocity fields (LES) in a vertical cross-section of the 10th cycle.

Fig. 11

 11 Fig. 11 Evolution of experimental ( ) and LES ( ) phase-averaged cylinder pressure during the whole cycle.

Figure 12 Fig. 12

 1212 Fig.[START_REF] Forsythe | Validation of a fluid-structure interaction model for a bileaflet mechanical heart valve[END_REF] Computed SGS viscosity normalized by fluid viscosity in the cylinder: volume-averaged ratio ( ), maximum ratio ( ).

Figures 13 and 14 Fig. 13

 1413 Figures 13 and 14 compare the measured data to the computed phase-averaged axial and vertical velocity fields in the median plane of the cylinder (plane 5 in Fig.1(b)) at four distinct times. The axial velocity fields are more representative of the flow motion induced by the valves whereas the vertical component is more representative of the movement of the piston. Note that the window of PIV measurement (Fig.1(b)) has been reduced to exclude zones where laser reflections or leaking oil issues were encountered in the near wall region. In the remaining window, the LES mean axial velocity component is in good agreement with the experimental data. The global shape and the level of velocities measured by PIV

Fig. 14

 14 Fig. 14 Comparison of LES (left) and PIV (right) results for the phase-averaged vertical velocity component (measurement noted 5 in Fig. 1(b)) at four instants.

Fig. 15

 15 Fig.[START_REF] Goryntsev | Large eddy simulation based analysis of the effects of cycle-to-cycle variations on air-fuel mixing in realistic DISI IC-engines[END_REF] Maximum in-cylinder pressure of cycle i + 1 (normalized by the maximum in-cylinder pressure of the mean cycle) versus the maximum in-cylinder pressure of the preceding cycle i; comparison between experimental (•) and LES (+) results.

Fig. 16

 16 Fig.[START_REF] Gourdain | High performance parallel computing of flows in complex geometries -part 1: methods[END_REF] The dashed line in (a) represents the line along which all of the profiles are plotted. (b)-(e) Comparison of axial velocity component profiles between the 25 instantaneous LES cycles ( ), the mean LES ( ), the mean experiment ( ) and the 95% confidence envelope of the 100 instantaneous experimental cycles ( ).

10 Fig. 17

 1017 Fig.[START_REF] Gourdain | High performance parallel computing of flows in complex geometries -part 2: applications[END_REF] Energy distribution of the first 10 most energetic POD modes.

  Fig. 18 Iso-surfaces of eigenfunctions of the POD first (black) and second (white) most energetic modes.

Table 1

 1 Characteristics of the SGEmac engine test bench. Crank Angle Degrees (CAD) are relative to compression Top-Dead-Center (TDC). The CAD of valve openings and closures correspond to experimental residual valve lifts of 0.1 mm.

		Units	Values
	Compression ratio	[-]	9.9
	Engine speed	[rpm]	1200
	Bore	[mm]	82
	Stroke	[mm]	83.5
	Connecting rod length	[mm]	144
	Intake Valve Opening (IVO)	[CAD]	350
	Intake Valve Closure (IVC)	[CAD]	-120
	Exhaust Valve Opening (EVO) [CAD]

Table 2

 2 Longitudinal mode frequencies predicted by the Helmholtz solver and by the LES for the intake line.

	Mode number name	location	Frequency [Hz]
				Helmholtz solver LES
	1	1/4 wave full intake line	108	110
	2	1/4 wave from valves to manifold	234	-
	3	3/4 wave full intake line	251	250
	4	5/4 wave full intake line	405	430

Table 3

 3 Convergence analysis on the number of experimental cycles used for ensemble averaging: differences with respect to reference quantities averaged over 100 cycles.

	CAD before		Mean velocity			RMS of velocity
	TDC	ε 10 ε 25 ε 50 ε 75	ε 10	ε 25	ε 50	ε 75
		%	%	%	%	%	%	%	%
	280	9.7 5.5 3.9	3	22.4 17.7 11.4 5.2
	240	8.4 6.2 2.3 1.7	16.4	9.3	5.2	1.9
	180	15	6.4 4.7 2.3	19.9	9.7	7.1	4.4
	100	17	5.6 3.6 1.3	20	12	7	2.3

Table 4

 4 Convergence analysis on the number of LES cycles used for ensemble averaging: differences with respect to reference quantities averaged over 25 cycles.

	CAD before		Mean velocity			RMS of velocity
	TDC	ε 5	ε 10	ε 15 ε 20	ε 5	ε 10	ε 15	ε 20
		%	%	%	%	%	%	%	%
	280	23.4 13.8 6.5 3.5	47.4 24.2 13.2 6.9
	240	17.5	7.4	5.2 2.2	38.8 20.5 11.4 5.3
	180	13.3	6.1	5	4.3	49.1 17.5	8.5	8.1
	100	9.6	6.4	5.1 3.7	58.4 26.1 12.7 6.5

Acknowledgements

The authors gratefully acknowledge the financial support provided by PSA Peugeot Citroën and Renault SAS in the framework of thesis researches. This work was granted access to the HPC resources of CINES (Centre Informatique National de lEnseignement Supérieur) under the allocation 2009-c2009026055 made by GENCI (Grand Equipement National de Calcul Intensif). The authors also acknowledge the financial support by ANR (Agence Nationale de la Recherche) under grant N.ANR-06-PDIT-007-01 SGEmac.


Aznam Yacoub
email: aznam.yacoub@lsis.org

Maamar El

Amine Hamri
email: amine.hamri@lsis.org

Claudia Frydman
email: claudia.frydman@lsis.org

Chungman Seo
email: cseo@rtsync.com

Bernard P Zeigler
email: zeigler@rtsync.com

To- Wards

Maamar Yacoub Aznam

TOWARDS AN EXTENSION OF PROMELA FOR THE MODELING, SIMULATION AND VERIFICATION OF DISCRETE-EVENT SYSTEMS

Keywords: DEv-PROMELA, Simulation, Formal Verification, Verification and Validation

d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Process Meta Language (PROMELA) is a well-known formalism proposed by [START_REF] Holzmann | Design and Validation of Computer Protocols[END_REF][START_REF] Holzmann | The SPIN Model Checker: Primer and Reference Manual[END_REF] for the modeling and the verification of concurrent systems by model-checking. Model-Checking [START_REF] Clarke | Design and Synthesis of Synchronization Skeletons Using Branching-Time Temporal Logic[END_REF]Emerson 1982, Baier and[START_REF] Baier | Principles of Model Checking[END_REF], and more generally formal verification techniques [START_REF] Huth | Logic in Computer Sciences: Modelling and Reasonning about Systems[END_REF], represent promising methods of verification. Because they are potentially able to explore all the states of a model, these methods appear as comprehensive and efficient methods. A verification model is focusing on the conceptual aspects of a design that are relevant to the properties one wants to verify. The meaning of that is that a verification model must be as simpler as possible, and PROMELA is designed to encourage the user to make abstraction of the computational aspects of the system under study. This also guarantees the efficiency of the model-checking algorithms. But, by enforcing these restrictions, the expressiveness of PROMELA is reduced. For example, the interactions between the components of a system can strongly depends on the values of data or on their timely coordination. Making strong assumptions about that or a too high-level abstraction can lead to the development of a non-efficient model.

In the opposite, Discrete-Event system Specification (DEVS) formalism [START_REF] Zeigler | Theory of Modeling and Simulation[END_REF]) is well-suited for lowlevel modeling and analysis of real systems. By providing a clean operational semantics, the DEVS formalism allows a clean interpretation of the model elements in the real world. The few restrictions enforced by DEVS increase its expressiveness. The other side of the coin is that it becomes impossible to apply formal proofs on these models. The challenge to introduce the DEVS semantics into a formalism of formal verification like PROMELA then seems to provide some benefits. First, it allows accurate modeling of discrete-event systems (DES) by using a formalism which supports formal verification. In other words, the purpose of such an approach is increasing the expressiveness of PROMELA without breaking its formal verification capabilities. Second, specifying a DES model into a syntactic formalism like PROMELA can make easier the translation from the conceptual model to a computerized simulation model. Indeed, transformation rules can easily be defined and verified between two syntactic formalisms. Third, the verification and validation of the conceptual DES models can be done through two combined methods, formal verification and discreteevent simulation, without needing to use multiple formalisms. The goal of this paper is to present preliminary work about an extension of PROMELA for modeling of DES. The section 2 is a quick overview of PROMELA. The section 3 presents our extension. Two example is given in section 4 to illustrate the verification process using the MS4Me environment [START_REF] Seo | DEVS Modeling and Simulation Methodology with MS4 Me Software[END_REF]) and the SPIN (Simple PROMELA Interpreter) model-checker [START_REF] Holzmann | Design and Validation of Computer Protocols[END_REF].

QUICK OVERVIEW OF PROMELA

PROMELA is a specification language with a semantics of executability. Its syntax is influenced by the Dijkstra's Guarded Command Language [START_REF] Dijkstra | Guarded commands, nondeterminacy and formal derivation of programs[END_REF]) and the C programming language defined by [START_REF] Kernighan | The C Programming Language[END_REF]. This fact makes that its use is relatively easy compared with other formal methods because its syntax is really close to any implementation languages. We are focusing in this section only on concepts interesting to the scope of this paper. The section 2.1 is about how specifying components of a system in PROMELA. The section 2.2 expose the semantics of PROMELA. We don't deliberately talk about the specification of the properties that one wants to verify on the model in this paper, because we are focusing on how increasing the accuracy of conceptual models.

System Specifications

A PROMELA system relies on three main types of objects: processes, data objects and messages. The components of the system are modeled by a finite set of instances of processes. The processes can communicate with each other thanks to different mechanisms: buffered messages, shared global variables or rendez-vous handshakes. Each process is a finite set of guarded commands called instructions. At any time t, only one instruction is executed without any assumptions about its duration. Note that a set of instructions can be labeled as atomic: in this case, these instructions are considered as a unique instruction, like any common atomic operations. Processes can be also prioritized, meaning that a process with a higher priority will always execute its instructions before other processes. Syntactically, a process is defined by a proctype block of instructions, as given in Program 1.

Instructions are divided into two categories: statements that modify the state of the system on the one hand, and control-flow instructions on the other hand. [START_REF] Natarajan | Outline for an Operational Semantics of PROMELA. The SPIN Verification System[END_REF]Holzmann 1997, Holzmann 2004). We can study it at two levels: on the one hand, at the process level, and on the other hand, at the whole program level.

Executability of a Process

A PROMELA process with a set L of statements is a finite state machine P = (Q, T, q 0 , F) where

 = { qi = (i, l1, …, lm, g1, …, gn, c1, …, co) ∈ ℕ x ∏ x ∏
x ∏ } is the finite set of states. Each state is defined by its id, the value of each local and global variable, and the value of each channel;  ⊂ x L x is the set of labeled transitions;  is the initial state;  is the set of final states.

Denote (q i ,q j ) ∈ Q 2 and l ∈ L (l is an instruction in the process). Then, t = (q i ,l,q j ) ∈ T iff the process can change its state from q i to q j by only executing l. This means that an instruction syntactically denotes two consecutive states. A transition can be executed only if it is enabled:

 l is a non-blocking instruction: an assignment, a conditional instruction with a satisfied guard, or any control-flow atomic instruction (else, skip, break, etc);  l is an asynchronous message sent over a nonfull channel;  l is an asynchronous message received from a non-empty channel;  l is an unblocking rendez-vous message.

If a transition is enabled and executed, the process changes its state from the source state to the target state by applying an appropriate action function. This action function modifies the values of the variables, the content of the channels, etc. If more of one instruction is enabled, the semantics engine randomly selects one of them and executes it.

Executability of a Program

A PROMELA program can be simply seen as the asynchronous product of the automata of each process. Then, a PROMELA program is defined by:

 ⊆ ∏
where is the set of states of the process P i ;  is the set of transitions; t = (si,l,sj) ∈ T if it exists a transition that changes the state of any of the processes by applying l;  is the initial state;

A transition t is enabled if its equivalent t i is enabled in the process pi. If two or more transitions are enabled, the semantics engine selects one of them at random and executes it (unless there is a defined priority between processes).

DISCRETE-EVENT PROMELA

The previous section shortly described the PROMELA semantics. As seen, PROMELA does not define either explicit time representation or events. But as stated by [START_REF] Tripakis | Extending PROMELA and SPIN for Real-Time[END_REF], the systems which can be modeled with PROMELA can be characterized as real-time systems.

The representation of time becomes thus important to develop a more efficient design. Many timed extensions of PROMELA [START_REF] Tripakis | Extending PROMELA and SPIN for Real-Time[END_REF][START_REF] Bosnacki | Discrete-Time PROMELA and SPIN[END_REF][START_REF] Nabialek | Translation of Timed PROMELA to Timed Automata with Discrete Data[END_REF]) and the algorithms for their verification were studied in the literature. In opposite of these approaches which try to deal with the algorithms of formal verification and with the PROMELA semantics, we propose to integrate the DEVS semantics into PROMELA. This approach has two main advantages. On the one hand, it gives to PROMELA the capability of modeling DES with a lower level of abstraction. The model is more accurate and formal verification can be applied only on interesting paths. Formal verification can also be used to check structural properties on the model. Moreover, while formal verification deals with finite models, discreteevent simulation can be used for simulation-based verification of behavioral properties on infinite models, validation and analysis. On the other hand, modeling DES with a syntactic language close to an implementation language make easier the translation from the conceptual model to the computerized simulation model. The sections 3.1 and 3.2 will recall concepts around discrete-event systems and the DEVS operational semantics. Then, the section 3.3 will introduce our new extension called Discrete-Event PROMELA (DEv-PROMELA).

3.1. Discrete-Event Concepts DES are a specific class of timed systems [START_REF] Zeigler | Theory of Modeling and Simulation[END_REF][START_REF] Zeigler | Multifacetted Modelling and Discrete-Event Simulation[END_REF]. A DES can be seen as an extension of a Moore Machine, in the sense that the outputs of a DES can depend only on its current state. [START_REF] Zeigler | Theory of Modeling and Simulation[END_REF] introduced two important concepts: the association of a lifespan to each state and the hierarchical composition. A DES evolves along the events that it emits or consumes, the distribution of events can be non-linear and time can be represented (through the lifespan concept) by any value of the ℝ -space. A DES thus relies on the following notions:

1. Each state has a lifespan whose the value is a real. When the lifetime of state expired, the current state changes according to the transition table. An output is then emitted; 2. When an input is consumed by a DES, its current state changes according to the transition table, regardless the current lifetime of the current state; 3. Each event are well-dated. If e1 and e2 are two events, thus e1 and e2 can always be ordered; 4. If two events overcome at the same time, then either they are probably equivalent events (e1 = e2) or they are prioritized; the priority between two events is well-defined (non-deterministic behavior are not permitted); 5. The state trajectory is stable between two successful events; 6. As a result of the previous points, transitions are characterized by their nature: internal transitions correspond to autonomous behaviors while external transitions allow the system to react to any external events.

Overview of the DEVS Formalism

A DEVS system is built upon DEVS atomic models which can be coupled to get a DEVS coupled model. A DEVS atomic model is then a small simulation unit defined by = ( , , , , , , ) where  X is the set of inputs;  Y is the set of outputs;  S is the set of sequential states;  : x → is the external transition function;

 : → is the internal transition function;  ∶ → is the output function;  ∶ → ℝ is the time advance function;  = {( , ) | ∈ , 0 ≤ dt ≤ ta(s)} is the set of total states.
Taking a DEVS atomic model A in a current state q. If the lifespan of q is ended (i.e. q = (s,ta(s))), the model changes its state to s'= (s). It emits the output y = λ(s). If A receives an event x, it change its state according to s'= (q,x). Otherwise, it stays in q, letting the time progress. Clocks can be used to define exactly the date of next events. The purpose of clocks is only to determine the next event by taking the minimum value of the dates of all next events. Missing concepts of events, state lifespan and type of transitions as viewed in the previous section are then introduced. Syntactically, a state is defined between two instructions in PROMELA. This means we can easily define the lifespan associated to each state and the type (internal or external) of the transitions associated to each instruction. We extend each statement with the following grammar described in the Backus-Naur Form: Each statement is prefixed by an event descriptor (in brackets) which defines the lifespan of the source state of the considering transitions and how transitions are triggered. For instance, look at these instructions:

[clt:3.0->emit:newa] a = a + 3;

(1)

[evt:newc] d = a + 9;

(2)

[clt:lifespan(a)->emit:newa | evt:newb ] b = a * 7;

(3)

(1) defines a state with a lifetime equal to 3 units of time, and an internal transition will emit the event "newa". By executing this transition, the assignment is executed.

(2) defines a passive state (with an infinite lifetime) and an external transition which will be enabled by receiving the event "newc".

(3) defines a state with a lifetime equal to "lifespan(a)" and two transitions: an internal one which will emit the event "newa" and an external one which will be enabled by consuming the event "newb". This third example shows how event descriptors can define transitions either as internal (with the clt command), or as external (with the evt command), or as a combination of multiple external events with probably one internal transition. A PROMELA statement can now label several transitions. Note that the default lifespan of each state is equal to infinity if an evt descriptor is defined and if there is no clt command. The clt command is optional and needed only to denote internal transitions. For convenient purposes, if a statement is not prefixed by an event descriptor, that is interpreted like it is prefixed by [clt:0 ->emit:silent]. The statement is executed without any delay by emitting the silent event. Priority descriptors on processes are also added. Indeed, the concept of DEv-PROMELA priority allows the ordering the events coming (resp. generating) from (resp. by) each process. This is a different meaning than the PROMELA initial concept of priority.

Semantics of a DEv-PROMELA Process

A DEv-PROMELA process P with a set L of statements is an automaton = ( , , , , , ) where

 = { si = ( , i, l1, …, lm, g1, …, gn, c1, …, co) ∈ ℕ x ∏ x ∏ x ∏ } is the finite set of states; 
is the finite set of events; E contains at least the silent event denoted ε;  ∶ → x E is the partial function that defines transitions labeled by a statement l ∈ L; note that it is a partial function because no internal transition can be defined for passive state.



∶ x E → is the partial function that defines external transitions labeled by a statement l ∈ L; note that it is a partial function because for a state q, there is not necessarily an external transition defined for each event e;  is the initial state;  is the set of final states.

We also define:

 : → ℝ ( ) ↦
is the lifespan function; Take a look at the Program 2. In the meaning of PROMELA, both guards are satisfied, meaning one of them is non-deterministically chosen. From another point of view, this means the model checker will verify the both paths, and this is the purpose of a good low-level verification model. But, if we take a look at the state space (Figure 3), we can clearly note that there is three cases. Two of them are unambiguous; concerning the intersection, the behavior is not well-defined.

Transposed to an event system, this model would describe two situations:

1. Denotes e 1 and e 2 two different events that overcome at the same date. The verification model does not care about which event is consumed at first; the two orders will be checked; 2. The lifespan of the state (x=2, y=2) is ended, but two behaviors are possible. In this case, the verification model is not well-defined and effort are put into the verification of paths that are maybe meaningless for the real system.

We mean that the non-determinism involved by the PROMELA control-flow structures doesn't come necessarily from the non-deterministic nature of the considered program, or from a high-level abstraction, but can be involved by incomplete specifications. This program can be interpreted in three manner:

1. if t1 ≠ t2, there is not ambiguity: the state concerned by the branches are not the same. Indeed, if t 1 is lesser than t 2 , T1 will be always triggered in first. If t2 is lesser than t1, T2 will be always triggered in first; 2. if t1 = t2, this model is non-deterministic.

Determinism is then enforced by always executing the first enabled branch T1 in the case (x=2, y=2).

This property is applicable to any branching structures.

Semantics of a DEv-PROMELA Program

A DEv-PROMELA program is a product of the automata of each process that composed the program, as in PROMELA. The difference is that, while in PROMELA the sequences of the states are interleaved (i.e. the final graph is composed by all the possible permutations of the statements), DEv-PROMELA takes into account only the execution paths given by the order of the events. Denote e 1 and e 2 , two events of the DEv-PROMELA system composed by two processes P1 and P2. If date(e 1 ) < date(e2), then the statement associated to e1 will be executed in first. Both processes update their state according to the action function and computes their local next event e 1 ' and e 2 '. The next event date is then given by min(e1' ; e2').

If date(e1) = date(e2), priority as defined in the previous section determines the next event between both events. Non-determinism is implicitly resolved for a precise instance: this does not mean that non-determinism doesn't exist in the global system. Because order of events can depend on the execution context, multiple branches could exist for a given global state. A verification algorithm would then have to check any possible orders of events. But in opposite of PROMELA, for a given order of events, one and only one path can be walked through.

It is also important to note that if the value of a global variable is changed, an implicit event is emitted. Other processes change their respective state. Their respective clock is however preserved, meaning their respective next event remains unchanged.

Fischer's Mutual Exclusion Protocol

The Fischer's Mutual Exclusion Protocol [START_REF] Abadi | An old-fashioned Recipe for Real-Time[END_REF]) looks like a very simple algorithm (Algorithm 4) for handling mutual exclusion: firstly, the process p checks whether another process either is already or wants to enter the critical section (line 4). If it is the case, the process stays in an active wait. Then, the process p declares his willing to enter the critical section before entering a sleep mode. When it wakes up, it checks whether another process is entered the critical section during its sleeping. If it is the case, the protocol restarts from the beginning, else the process can enter the critical section. Fisher's protocol can be seen as a timed system, and also a timed event system.

Using DEv-PROMELA, it is easy to obtain a verification and a simulation model (Program 6) of this algorithm. Because of the system can be considered as closed, almost of the transitions correspond to an autonomous behaviors which depend only on the current state. Priority ensures that each process will have the hand in turn. When a process is sleeping, it is blocked until the delay is exactly expired or until another process changes the value of the id variable. This is modeled in l.16: the model stays in this state during delay units of time or until it receives the changepid event. In other words, the next event for a process will be at getCurrentDate+delay after the execution of line 14.

If both processes wake up at the same time, priority will decide of the next event to proceed. Remember that the execution of this DEv-PROMELA program is deterministic.

The structural preservation property of DEv-PROMELA is out-of-scope of this paper. This property allows the generation of an untimed PROMELA verification model which has the same structural properties than the DEv-PROMELA model (Program 5). Thanks to that model, two properties can be easily verified:

1. "Two processes cannot be in critical section in same time." (□ ! ( ∧ ))" 2. "Is there cases in which two processes cannot be in critical section in same time ?" (◊ ! ( ∧ ))

While the property 2 is well-verified by the model, the SPIN model-checker returns a violation error as expected for the property 1. Indeed, a process P can execute the l.11 and l.14, then a process Q can execute the l.11 and l.14 and enters the critical section while P is always in it. This is well-known that the timed constraint delay > C (where C is the longest time that a process may take to perform a step while trying to enter the critical section) must be met to guarantee the mutual exclusion by this protocol. This only result could enforce the designers to design solutions to prevent such an error. But if the system is well specified, for example if the duration taken to execute each instruction is quantified, a simulation-based verification can bring another interpretation. Remember that the DEv-PROMELA semantics is the same than the DEVS semantics, meaning a DEv-PROMELA model is a DEVS model. As a consequent, a DEv-PROMELA model is also a simulation model. We encode the Program 6 into a DNL format given in Appendix A by translation; then, we simulate it into the MS4 Me Environment (Figure 2). A simulation-based verification is then performed. In the context in which two processes exactly acts in the same manner and with the semantics of Classic DEVS, the analysis reveals that whatever the value of the delay, two processes can't be in the critical section at the same time.

Alternating Bit Protocol

We also use DEv-PROMELA to model the Alternating Bit Protocol (APB) [START_REF] Bartlett | A note on reliable full-duplex transmission over half-duplex lines[END_REF], which has been used to illustrate the analysis capabilities of PROMELA [START_REF] Holzmann | The SPIN Model Checker: Primer and Reference Manual[END_REF]).

CONCLUSION AND FUTURE WORK

We described in this paper some preliminary work to develop a new extension of PROMELA for the modeling and simulation of DES. This is done by introducing the DEVS operational semantics into PROMELA, without breaking the verification capabilities of PROMELA. A DEv-PROMELA model is then a verification model and a simulation model which can be verified by combining formal verification for structural properties and simulation-based verification for behavioral properties.

The other advantage is that we are providing a syntactic formalism for DEVS modeling. This can help to implement computerized models from conceptual models. This maybe also enables the formal verification of the simulation models by using rewriting rules (for example from DEv-PROMELA to DNL).

The major drawbacks of our methodology is that DEv-PROMELA is limited to the Classic DEVS and its subclasses. Parallel DEVS cannot be modelled with this approach due to the limitations of PROMELA. But it can be interested to study how introducing the semantics of Parallel DEVS into PROMELA. Future work also concerns providing a formal proof that a DEv-PROMELA model is well a DEVS model, and working on limiting the formal verification algorithms only to the paths which are really expressed by the DEv-PROMELA models, without generating all the state space of the PROMELA equivalent.

  <proctype decl> ::= "[" priority "=" <int> "]" <proctype> <event stmnt> ::= "[" <timed trans> "]" <stmnt> | <stmnt> <timed trans> ::= <clt expr> | <evt expr> | <clt expr> "|" <evt expr> <clt expr> ::= "clt:" <real expr> "->emit:" <evt val> <evt expr> ::= "evt:" <evt val> [ <op> <evt expr> ] <op> ::= " |" <evt val> ::= <mtype> | "silent" <real expr> ::= <real> | "infinity" | /* Any C-function returning a real value */



  = { = ( , ), ∈ and 0 ≤ dt ≤ ta(s) } is the set of total states;  = { = ( , 0), ∈ }  = { = ( , ( )), ∈ } Take s a state, and l a statement with an event descriptor.If l denotes an internal transition, l is enabled and (s) is triggered only if the value of the current clock of the process is equal to the lifetime ta(s). The event associated to the transition is emitted and the next event date for the process is computed by: de = getCurrentDate + ta(s') if (s'; e) = (s).If l denotes an external transition on an event e, l is enabled if the process receives the event e. If the next event of the current process is at the same date than another external event, the priority statement resolves the conflict. If ta(s) = ∞, all of the statements related to s are blocked. Only external transitions can be enabled in this case. Because there is at most one internal transition per state and because there is a priority statement, a DEv-PROMELA process is deterministic. What about the non-determinism expressed in the PROMELA control-flow structures ?

Figure 1 :

 1 Figure 1 : State space of the Algorithm 2. T1 denotes the transition labelled by the statement line 4, and T2 the transition labelled by the statement line 5.

  extension to model two examples of systems. First is the well-known Fischer's Mutual Exclusion Protocol. Second concerns the Alternate Bit Protocol with Finite Queue and Infinite Queue. Verification of properties were done through combined formal verification using the classic SPIN modelchecker, and simulation-based verification using the MS4Me Environment. Because the DEv-PROMELA semantics is the same than the one of DEVS, the translation from DEv-PROMELA specifications to DNL format were easy.

  Property 1 is thus verified by this DEv-PROMELA model. How interpreting the difference between the results of formal verification and simulation ? On the one hand, the formal verification gives first information about the correctness of the model, on a high-level abstraction. For example, the formal verification ensures the absence of deadlocks. It can help designer to fix the design before considering other forms of testing. On the other hand, the simulation-based verification gives more precise information to the designer to understand what really happened. Because the DEv-PROMELA model is at lower level of abstraction than the PROMELA model, the designers can decide that the error found by the model checker is outside the scope of the model or not. Another advantage is that discrete-event simulation can be also used for validation at the same time. Designer can then get information about the validity of this model against its intended use.

  The only things we have to do is to model the time T to transmit a bit, and the time T1 up to that we can consider a message is lost. The obtained model does not really differ from the PROMELA one. But what happened if we want take into account a queue which stores messages before sending them ? As shown by[START_REF] Zeigler | Combining DEVS and Model-Checking: Using System Morphisms for Integrating Simulation and Analysis in Model Engineering[END_REF], the capacity of the queue has an impact on the behavior of the system. Modeling the queue in PROMELA certainly increases the size of the state space, and reduces the efficiency of the formal verification algorithms. But DEv-PROMELA can be used in two manner: as a verification model, assumptions can be verified in idealized conditions (in this case, by abstracting the capacity of the queue or by considering there is no queue), and as a simulation model. Then, the DEv-PROMELA model of the APB can be coupled with any simulation model of queues. The simulation-based verification then extends the formal verification by exploring scenarios outside the scope of the formal verification.

Figure 2 :

 2 Figure 2 : State Diagram of the DEv-PROMELA model of Fischer's Algorithm in MS4 Me Environment.

  Assignments involve local and global variables, whereas communication statements involve global buffered channels. Control flow instructions are classical conditional and loop structures. They allow the selection of the next statement among different branches regarding a guard. Because PROMELA processes are nondeterministic, if several guards are satisfied, the next instruction is randomly selected. If none of them is satisfied, the control flow structure is blocked.PROMELA provides two special guards called else and timeout: if these guards are present in a control flow structure, and if this one is blocked, these instructions are then executed. Data in PROMELA is represented by local and shared variables. Local variables are relative to the process only in which they are declared, whereas global variables are shared by all the processes. A variable is characterized by its value and its type, among either all the PROMELA scalar basic datatypes as given in Table1, or any finite combinations (structures) or finite arrays of these types.

	Table 1 : Basic PROMELA Datatypes
	Type	Size (bit) Range of Values
	bit, bool	1	[0;1]
	byte	8	[0;255]
	mtype (constants)	8	[0;255]
	short	16	[-2 16 ;2 16 -1]
	int	32	[-2 32 ;2 32 -1]
	2.2. PROMELA Semantics
	PROMELA centers on a semantics of executability

   Y is the set of outputs;  D is the name set of components;  M is the set of atomic components;  EOC is the set of external output coupling;  EIC is the set of internal input coupling;  IC is the set of internal coupling;  Select is the tie-breaking function to choose the next event from the set of simultaneous events.Informally, a DEVS coupled model is composed by atomic models. The date of the next event is computed by getting the minimum value among all the dates of all the next events. By processing an event, each component is updated according to its transition functions. If multiple events occur at the same date, the Select function chooses the next event to proceed.

	3.3. Building Discrete-Event PROMELA
	Introducing a new semantics in PROMELA needs to
	define new syntactic elements to model the concepts we
	are adding. The section 3.3.1 introduces these new
	syntactic elements. The section 3.3.2 and 3.3.3 define the
	meaning a DEv-PROMELA model.
	A DEVS coupled model is defined by	=
	( , , , { },	,	, ,	) where
	 X is the set of inputs;

3.3.1. Syntactic Rules

A new datatype called real is introduced. The purpose of real variables is to represent infinite and unbounded real values. Real variables can be local or global, and their declaration does not differ from any scalar variables: real i,j,k; Real variables can be combined in complex structures or in arrays, without any restriction. An implicit local clock is associated to each process. Consequently, like in Timed PROMELA, introducing a specific type for clocks is not needed. Each process is ensured to have a clock. The clock valuation can be accessed through a specific getCurrentDate instruction.

Proceedings of the European Modeling and Simulation Symposium, 2015 978-88-97999-57-7; Affenzeller, Bruzzone, Jiménez, Longo, Merkuryev, Zhang Eds.

ACKNOWLEDGMENTS

This work is a part of the R&D project "MAGE", from French "Investing for the Future" national program.

APPENDIX A: DEV-PROMELA MODEL OF FISCHER'S ALGORITHM ENCODING IN THE DNL FORMAT

use delay with type double and default "0.1"! use id with type int and default "1"! use pid with type int and default "0"! generates output on silent! generates output on changedpid1! generates output on changedpid0! accepts input on changedpid0! accepts input on changedpid2! to start hold in S1 for time 0.1! after S1 output silent! from S1 go to S2! when in S1 and receive changedpid0 go to S1! when in S1 and receive changedpid2 go to S1! external event for S1 with changedpid0 <%pid = 0;%>! external event for S1 with changedpid2 <%pid = 2;%>! internal event for S1 <%if(pid==0) holdIn("S2",0.1); else holdIn("S1",0.1);%>! hold in S2 for time 0.1! after S2 output changedpid1! from S2 go to S3! when in S2 and receive changedpid0 go to S2! when in S2 and receive changedpid2 go to S2! external event for S2 with changedpid0 <%pid = 0;%>! external event for S2 with changedpid2 <%pid = 2;%>! internal event for S2 <%pid = id;%>! hold in S3 for time delay! when in S3 and receive changedpid0 go to S1! when in S3 and receive changedpid2 go to S1! external event for S3 with changedpid0 <%pid = 0; holdIn("S1",0.1);%>! external event for S3 with changedpid2 <%pid = 2; holdIn("S1",0.1);%>! after S3 output silent! from S3 go to S4! internal event for S3 <%if(pid==id) holdIn("S4",0.1); else holdIn("S1",0.1);%>! hold in S4 for time 0.1! after S4 output changedpid0! from S4 go to S1! when in S4 and receive changedpid0 go to S4! when in S4 and receive changedpid2 go to S4! external event for S4 with changedpid0 <%pid = 0;%>! external event for S4 with changedpid2 <%pid = 2;%>! internal event for S4 <%pid = 0;%>!