Ionel Nistor

Olivier Pantalé
email: olivier.pantale@enit.fr

Serge Caperaa

Numerical propagation of dynamic cracks using X-FEM

Keywords: partition of unity, extended finite element method, dynamic crack propagation, cohesive model. MOTS-CLÉS : partition de l'unité, méthode des éléments finis étendue, propagation dynamique de fissures, modèle cohésif

This paper presents an application of the eXtended Finite Element Method for numerical modeling of the dynamic cracks propagation. The numerical cracks representation is adapted to the time-dependent mechanical formulation, using the Heaviside step function for completely cutted elements and the cohesive model for crack tips. In order to find the propagation parameters, a crack evolution model is proposed. The numerical implementation is achieved in new explicit FE module. A numerical example is proposed for proving the computational efficiency of this new module.

RÉSUMÉ. Cet article présente une application de la méthode des éléments finis étendue à la simulation numérique de la propagation des fissures dynamiques. La représentation numérique des fissures a été adaptée à la dépendance au temps de la formulation mécanique en utilisant la fonction Heaviside pour les éléments complètement coupés et le modèle cohésif pour la tête de la fissure. Un modèle d'évolution de la fissure est également proposé. L'implémentation numérique a été réalisée dans un nouveau module de calcul éléments-finis explicite. Un exemple numérique est présenté pour démontrer l'efficacité de ce nouveau module de calcul.

Introduction

The numerical modeling of the discontinuous structures under dynamic loads is a major challenge of the last decades. Several techniques were developed in order to predict the behavior of such discontinuities, cracks or shear bands.

One of the earlier technique used for discontinuities representation is the remeshing. Based on the classical FEM, the mesh of the analysed structures is re-build at each time step in order to propagate the crack over the elements edges/faces. Because of simplicity, this technique is implemented in several commercial FE codes and the later developments improved its features (the re-meshing algorithm acts only around the crack tips). However, several drawbacks remain: the mesh sensitivity on the direction of the crack propagation or the managing of the output data for the geometrical points around the cracks. Therefore the application of this technique to dynamic crack propagation remais quite difficult but some later contributions published on this subject proposed solutions for the instability problems. So, Réthoré et al. [START_REF] Réthoré | A stable numerical scheme for the finite element simulation of dynamic crack propagation with remeshing[END_REF] proposed a technique called "balance recovery method" which guarantees both numerical stability and accuracy for any type of projection used in dynamic simulations with varying meshes.

Later, in order to avoid the insertion of supplementary nodes during the re-meshing step, the inter-elements methods with cohesive laws appeared. So, cohesive segments are placed along the element interfaces as the crack advances. This approach was followed by [START_REF] Xu | Numerical Simulation of Fast Crack Growth in Brittle Solids[END_REF], [START_REF] Camacho | Computational modeling of impact damage in brittle materials[END_REF][START_REF] Ortiz | Finite-Deformation Irreversible Cohesive Elements for Threedimensional Crack-Propagation Analysis[END_REF]Pandolfi et al., 1999) and later by [START_REF] Remmers | Simulation of Fast Crack Growth Using Cohesive Segments[END_REF]. The efficiency of this technique has been proved for the simulation of dynamic crack propagation but the advancing direction is limited to the elements edges, and therefore, sensitive to the mesh orientation. However, the cohesive concept was developed for modeling the crack-tip area through some other methods, as one can see later in this paper.

Another class of fracture analysis techniques concerns the embedded discontinuity methods, inspired by the work of Ortiz, Leroy and Needleman [START_REF] Ortiz | A Finite Element Method for Localized Failure Analysis[END_REF]. The term "embedded discontinuity" was used by [START_REF] Belytschko | A Finite Element with Embedded Localisation Zones[END_REF]. This one considers that the discontinuity is a high-strain narrow band between two weak discontinuity lines at level element. The same approach was used for strong discontinuities by [START_REF] Dvorkin | Finite Elements with Displacements Interpolated Embedded Localization Lines Insensitive to Mesh Size and Distorsions[END_REF] and [START_REF] Simo | An Analysis of Strong Discontinuities Induced by Strain Softening in Rate-Independent Inelastic Solids[END_REF]. An extensive review of this class of methods is given by [START_REF] Jirasek | Comparative Study on Finite Elements with Embedded Discontinuities[END_REF].

The dependence of the crack advancing direction with mesh orientation remained one of the major drawback of the above mentioned methods foccusing the major part of the efforts in order to find a solution. The framework of the Partition of Unity Method (PUM), developed by [START_REF] Babuska | The Partition of Unity Method[END_REF] offered a solution. It states that the new function space, combination between regular shape functions and enriching ones, is able to represent the crack kinematics, without re-meshing. Using the features of PUM, the eXtended Finite Elements Method (XFEM) was developed as a method based on the enrichment of the standard displacement approximation with additional functions. Discontinuous partition of unity enrichment were firstly used by [START_REF] Belytschko | Elastic Crack Growth in Finite Elements with Minimal Remeshing[END_REF] who used the near-tip field enrichment and re-meshing technique for the rest of crack. The original term XFEM was introduced by [START_REF] Moës | A Finite Element Method for Crack Growth without Remeshing[END_REF] who developed the method proposing the introduction of a step function enrichment for the elements completely cutted by the crack.

This new approach was constantly improved and applied for different fracture problems: [START_REF] Daux | Arbitrary Branched and Intersecting Cracks with the Extended Finite Element Method[END_REF] for branching cracks, [START_REF] Sukumar | Extended Finite Element Method for Three-Dimensional Crack Modelling[END_REF] for threedimensional cracking, [START_REF] Moës | Extended Finite Element Method for Cohesive Crack growth[END_REF] for cohesive crack growth, [START_REF] Belytschko | Dynamic Crack Propagation Based on Loss of Hyperbolicity and a New Discontinuous Enrichment[END_REF][START_REF] Belytschko | Singular enrichment finite element method for elastodynamic crack propagation[END_REF] and [START_REF] Réthoré | An energy-conserving scheme for dynamic crack growth using the eXtended finite element method[END_REF] for elastodynamics crack propagation and, recently, [START_REF] Elguedj | Appropriate extended functions for X-FEM simulation of plastic fracture mechanics[END_REF] proposed the appropriate extended functions for modeling the plastic fracture using XFEM.

The objective of this work is to apply the XFEM approach for the simulation of the dynamic crack analysis and to implement it in a home-made FE code. This paper presents the mains steps of this work and is organized as follows. The crack representation fundamentals and the mechanical formulation according to the XFEM time-dependent model are described in Section 2. Next, a dynamic crack propagation criterion based on a cohesive zone concept is proposed. Some essential aspects concerning the numerical implementation are presented in Section 4 of the paper. Some numerical results for a dynamic crack analysis example are presented in Section 5 in order to illustrate the capabilities of the new code. Finally, some conclusions and prospects are reported in the last Section of this paper.

Crack representation in XFEM time-dependent model

As mentioned before, in order to represent the cracks in XFEM, the framework of PUM was used for the enrichment of the classical displacement-based finite element approximation. Therefore, the additional degrees of freedom are introduced for the nodes whose nodal shape function support are cutted by the crack. In the initial development of XFEM [START_REF] Moës | A Finite Element Method for Crack Growth without Remeshing[END_REF], for quasi-static analysis, two different types of enrichment are considered: the Heaviside step function for the elements completely cutted by the crack and the Westergaard asymptotic near-tip function for the elements containing the crack-tip.

The time-dependence of the dynamic analysis problem makes more difficult the implementation of this mixed enrichment in XFEM, therefore, another solution, more simple, was proposed recently [START_REF] Belytschko | Dynamic Crack Propagation Based on Loss of Hyperbolicity and a New Discontinuous Enrichment[END_REF]. In this approach the crack-tip position is restricted to the elements edges, i.e. the crack-tip advances from edge to edge, avoiding in this way the enrichment of the near-tip area with the Westergaard asymptotic functions. So, only the Heaviside step function is used for enriching the classical displacement-based finite element approximation of all elements being cutted by the crack. This is the approach that we also adopted in our work. The crack-tip areas have to be treated by a cohesive model, as presented in the next Section of the paper. The Heaviside step function H( -→ X ) allows to represent the discontinuity in the displacement field cutted by the crack. The Heaviside step function takes the value +1 above the crack and -1 below the crack (considering the crack advancing direction) 1:

H( -→ X ) = +1 if -→ X - -→ X * -→ e n > 0 -1 otherwise [1]
where -→ X is the sample point, -→ X * is the closest point to -→ X on the crack and -→ e n is the unit outward normal to the crack at -→ X * . The additional degrees of freedom introduced here are attached to the nodes belonging to cutted elements. Therefore, the new discontinuous displacement field -→ u h for a N nodes mesh, including N Γ enriched nodes, is approximated by:

-→ u h ( -→ X ) = I∈N ϕ I ( -→ X ) -→ u I + I∈NΓ ϕ I ( -→ X )H( -→ X ) -→ a I [2]
where -→ u I corresponds to the classical degrees of freedom, -→ a I to the enriched ones and ϕ I ( -→ X ) are the shape functions of the element.

Figure 1. Notations for the 2D discontinuous body

Further, the equations of elastodynamics are recalled here. A cracked homogeneous domain Ω in the current configuration is considered (as shown in Figure 1). A crack is represented by the curve Γ c , a traction t is applied on the Neumann boundary Γ t and ū is the applied displacement on the Dirichlet boundary

Γ u , Γ u ∪ Γ t = Γ, Γ u ∩ Γ t = 0.
The strong form of mechanical equations, in terms of the Cauchy stresses, for the current configuration, can be written as follows:

∂σ ij ∂x j + ρb i -ρ •• u i = 0 in Ω [3]    n i σ ij = t i on Γ t n j σ - ji = n j σ + ji = τ on Γ c u i = ūi on Γ u [4]
where ρ is the current density, b is the body force per unit mass, τ is the cohesive traction across the crack and σ is the Cauchy stress tensor. The weak form of the momentum equation is then:

Ω δu i ρ •• u i dΩ = Ω δu i ρb i dΩ + Γt δu i ti dΓ t - Γc δu i τ i dΓ [5] - Ω\Γc ∂ (δu i ) ∂x j σ ij dΩ
where u i is the trial displacement field (see [2]), δu i is the test displacement field and by δu i is denoted the displacement jump across the crack. This last quantity is depends only on the magnitude of the enriched degrees of freedom (see [19]).

The equilibrium discrete equations for dynamic analysis with XFEM are obtained from Equation [5] using the standard Bubnov-Galerkin procedure, substituting trial and test displacement fields and its derivatives:

M uu IJ M ua IJ M au IJ M aa IJ •• u J •• a J + F int iI Q int iI - F ext iI Q ext iI = 0 [6]
with

F int iI = Ω\Γc ∂ϕ I ∂x j σ ij dΩ F ext iI = Γt ϕ I ti dΓ + Ω ϕ I ρb i dΩ [7] Q int iI = Ω\Γc ∂ (ϕ I H) ∂x j σ ij dΩ [8] Q ext iI = Γt (ϕ I H) ti dΓ + Ω (ϕ I H)ρb i dΩ - Γc 2ϕ I τ i dΓ [9] M uu IJ = Ω ρϕ I ϕ J dΩ M ua IJ = Ω ρϕ I (ϕ J H)dΩ [10] M aa IJ = Ω ρ(ϕ I H)(ϕ J H)dΩ [11]
where in the global mass matrix the mixed and enriched terms are presented denoted by index ua, au and aa, respectively. For the numerical integration the consistent mass matrix will be used since the terms linked to the enrichment degrees of freedom obstruct its lumping. A discution about this choice is welcomed here. Several authors treating dynamic crack propagation met this problem: using the consistent mass matrix or the lumped one in case of the presence of enriched degrees of freedom. [START_REF] Belytschko | Dynamic Crack Propagation Based on Loss of Hyperbolicity and a New Discontinuous Enrichment[END_REF] retained the terms linked to the enriched degrees of freedom from consistent mass matrix because that was unable to determine an effective diagonal mass matrix for these ones. [START_REF] Borst | Mesh-independent discrete numerical representations of cohesive-zone models[END_REF] considering that lumping a such mass matrix leads to the lost of the essential information on the coupling of the regular and enriched degrees of freedom. For the numerical representation of cohesive-zone models this fact gives rise to the spurious transmission of stress waves through cohesive surfaces. Recently, [START_REF] Menouillard | Efficient explicit time stepping of the eXtended Finite Element Method (X-FEM)[END_REF] introduced a lumped mass matrix for enriched elements. The stability of this technique is proved showing that the critical time step of an explicit integration scheme does not tend to zero if the crack reaches the elements boundaries. Nevertheless, in the particulary case of Heaviside step function, the coupling terms are not taken into account with this lumping technique. Our choice, to use the consistent mass matrix are reasoned by the presence of a cohesive model and acceptable in terms of computational cost for relative small number of total degrees of freedom. For avoiding the diminution of the critical step time, we ignored the cracked elements in the numerical algorithm for searching its value (as presented in Section 4).

Dynamic crack propagation model

For a complete analyse of a cracked structure using XFEM one needs a dynamic crack propagation model, aside the discrete equations system given by [6]. This one must supply the answers to three essential questions: when the crack advances, in which direction and how fast?

Concerning the first question, in our development of XFEM for dynamic analysis, the choice of the crack propagation model is influenced by the enrichment type adopted. We use the cohesive area concept for modeling the crack tip and therefore the crack propagation model is based on this one. This approach was already used by [START_REF] Camacho | Computational modeling of impact damage in brittle materials[END_REF] for a lagrangian finite element method, [START_REF] Moës | Extended Finite Element Method for Cohesive Crack growth[END_REF] for quasi-static crack growth with XFEM and [START_REF] Belytschko | Dynamic Crack Propagation Based on Loss of Hyperbolicity and a New Discontinuous Enrichment[END_REF] for dynamic analysis with XFEM.

We called this model "cohesive damage model" and the main idea is to assume a cohesive zone located in front of the crack tip, in the elaboration zone of crack (Figure 2a) between the mathematical crack tip (i.e. the point where crack opening displacement vanishes) and the physical crack tip (i.e. the point of complete separation of crack lips). A cohesive law is considered along this zone. This one is linking two parameters of the model: the critical opening displacement corresponding to physical crack tip for a complete separation, noted δ crit , and the maximal value of cohesive stress at the mathematical crack-tip, noted σ max . The magnitude of the cohesive stress decreases as the crack opening displacement increases along the cohesive zone and the relationship between these two parameters, widely accepted in literature is the linear one (Figure 2b), the area under this curve matching the fracture energy:

G F = σ max δ crit 2 [12]
As one can see, considering this shape of the cohesive law, the cohesive crack model is completely defined by the fracture energy G F , the maximal cohesive traction values of σ crit and the cohesive zone length, Λ, given as material properties. The crack will propagate if the crack opening displacement exceeds the critical value, δ crit , inside of cohesive zone. In the next Section, the numerical evaluation of this parameter is detailed. This one gives the answer to the first question above.

The propagation direction is determined using the maximum circumferential stress criterion. The numerical searching algorithm used is presented in the next Section.

The particular approach adopted for this implementation of XFEM, restricting the position of the crack tip to the element edges, provides the crack propagation speed automatically since the crack tip is allowed to advance one element at a time.

Numerical implementation

Numerical implementation of XFEM for dynamic analysis was achieved using the framework of the DynELA explicit finite element code [START_REF] Pantalé | An Object-Oriented Programming of an Explicit Dynamics Code: Application to Impact Simulation[END_REF]. This FEM code is entirely programmed in C++ within an Object-Oriented Programming (OOP) approach. This feature presents a very well defined mechanism for modular design and re-use of code. In our case, this allows the development of the new module Dy-naCrack, intended for XFEM, without the need of re-implementing the major part the code. New classes have been created for the support of XFEM particularities and other have been specialized by the use of inheritance from DynELA classes. Further, we only present here some of the the particularly characteristics of this XFEM numerical implementation. A complete description of DynaCrack, is available in [START_REF] Nistor | Identification expérimentale et simulation numérique de l'endommagement en dynamique rapide: application aux structures aéronautiques[END_REF]. Currently, this module performs the elastodynamics analysis using the XFEM for bidimensional models presenting a single crack without branching. Only a 4-node full integrated quadrilateral element is implemented. Several improvements are planned in the next future, especially the extension to the propagation of dynamic cracks in plastic materials.

The first important challenge for the numerical implementation of XFEM is to achieve the numerical integration of the discrete equation over the domains of the cutted elements. The widely accepted solution in this case is the partitioning of these elements instead of re-meshing, since the re-meshing concept has some drawbacks as presented in the first Section of this paper. The partitioning procedure in XFEM is achieved for numerical integration purpose only, and no additional degrees of freedom are introduced into the discrete space. For 2-D models, the partitioning of cutted elements into a set of triangles was proposed by many authors [START_REF] Moës | A Finite Element Method for Crack Growth without Remeshing[END_REF][START_REF] Sukumar | Modelling Quasi-Static Crack Growth with the Extended Finite Element Method. Part I: Computer Implementation[END_REF][START_REF] Dolbow | Enrichment of Enhanced Assumed Strain Approximations for Representing Strong Discontinuities: Addresing Volumetric Incompressibility and the Discontinuous Patch test[END_REF], and several numerical integration options were presented also. Special (higher-order) quadrature rules have to be used on the elements that are partitioned in this way [START_REF] Moës | A Finite Element Method for Crack Growth without Remeshing[END_REF] like the use of a six-point integration rule. The solution that we adopted in our approach is to subdivide the two sub-domains separated by the crack for the cutted elements into sub-quadrilaterals, as illustrated in Figure 3. The main reason for this partitioning solution is related to the numerical integration accuracy. The numerical computation of the requested quantities as stiffness matrix on the sub-quadrilaterals is done using the same integration points as for an entire element of the mesh. In the same time, bilinear shape functions, used to support the quantities to integrate, are more appropriated to this Gauss integration rule than to a three or six-point integration rule on triangular domains.

Another important consequence of the used of XFEM description concerns the management of the "enriched" elements contribution to global matrices and vectors. We call enriched elements, the ones containing at least one enriched node. So, as one can see for the example in Figure 4, an XFEM mesh contains two types of enriched elements (beside the normal ones): completely cutted elements having all nodes enriched (in dark grey) and partially enriched elements supporting from one to three enriched nodes (in light grey).

Figure 4. XFEM enriched mesh with a crack

The size of the elementary matrices for these enriched elements are different depending on the number of degrees of freedom (classical and enriched) supported by the nodes. This makes the assembling procedure of the global matrices or vectors more difficult, taking into account the evolution of the elements status, too. As the crack advances through the mesh, some "normal" elements can became "enriched" ones. In order to achieve the assembling of the global matrices and vectors, the local and global positions of each degree of freedom in the linear algebraic system must be known at every time step. So, in DynaCrack, a special procedure was developed allowing a dynamic bi-directional mapping between local position of a degree of freedom and the global one. Up-dated at each integration time step, this procedure is useful for post-processing too.

For the numerical integration of the mechanical discrete equations, a displacements-based formulation was chosen. On behalf of the XFEM approach, for the enriched nodes, the displacements are composed of a continuous -→ u c (t) and an enriched -→ u e (t) contribution that takes into account the sign of the Heaviside function:

-→ u (t) = -→ u c (t) + -→ u e (t) [13]
An explicit processor similar to the one used by DynELA FEM code [START_REF] Pantalé | Parallelization of an Object-Oriented FEM Dynamics Code: Influence of the strategies on the SpeedUp[END_REF] was implemented using the Chung-Hulbert [START_REF] Hulbert | Explicit Time Integration for Structural Dynamics with Optimal Numerical Dissipation[END_REF] integration scheme:

•• -→ u n+1 = M -1 --→ F ext n - --→ F int n -α M •• -→ u n 1 -α M [14] • -→ u n+1 = • -→ u n +∆t n+1 (1 -γ) •• -→ u n +γ •• -→ u n+1 [15] -→ u n+1 = -→ u n + ∆t n+1 • - → u n +∆t 2 n+1 1 2 -β •• -→ u n +β •• -→ u n+1 [16]
The main feature of this explicit integration algorithm is the presence of some numerical dissipation defined by the values of the characteristic parameters α M , β and γ. The values of these parameters are given by the following relations:

α M = 2ρ b -1 1 + ρ b β = 5 -3ρ b (1 + ρ b ) 2 (2 -ρ b ) γ = 3 2 -α M [17]
where ρ b ∈ [0, 1], the spectral radius, defines the dissipation of the integration algorithm: ρ b = 1 leads to a conservative algorithm, since the lower value leads to a damping one. In the applications presented here, the conservative algorithm is considered. When the crack advances for t n+1 the new enriched degrees of freedom are added to the model and initialized to zero. This scheme for according the state vectors size, simillary to the one introduced by Réthoré et al. [START_REF] Réthoré | An energy-conserving scheme for dynamic crack growth using the eXtended finite element method[END_REF], has proved its stability.

For computing the time-step ∆t two approches were studied: the first one related to the maximal modal frequency of the structure and the second one related to the length of the smallest element of the domain and to elastic properties of the material. A lower computational cost was obtained with the second solution and finally the following relation was used:

∆t = f l e c d [18]
where f is a factor taken 0.82, l e is the caracteristic length of the smallest element of the domain and c d is the dilatational speed of the material. Note that the algorithm research for l e not take into account the elements cutted by the crack in order to avoid the diminution of the time step when the crack passes closely by enriched nodes.

Specific numerical algorithms were developed for implementing the cohesivedamage crack propagation model in DynaCrack. The evaluation of the advancing criteria, illustrated in Figure 5a, is a quite easy task since the discontinuity δ i in the new-introduced "numerical crack-tip" C n is evaluated for the last cutted element (E b in the figure), by:

δ i = 2 4 J=1 ϕ J (C n ) a J [19]
where a J are the enriched degrees of freedom. If the discontinuity jump given by [19], which represents the crack opening displacement value at the numerical crack tip, exceeds the critical value δ crit , the crack is propagated over the next element. As the numerical crack tip is always located on an element edge, the displacement is fixed, therefore, the time needed for δ i ≥ δ crit determines the advancing speed.

The crack advancing direction is given by the searching algorithm shown in Figure 5b. The maximum circumferential stress criterion sets that a crack should propagate from its tip in the direction θ c of maximal circumferential stress σ θθ . The evaluation of σ θθ is performed on a half-circle centered on the crack tip. Erratic crack paths due to the use of this criterion were reported [START_REF] Belytschko | Dynamic Crack Propagation Based on Loss of Hyperbolicity and a New Discontinuous Enrichment[END_REF] proved by our numerical simulations too. With the use of an appropriate value for the searching radius, validated for fracture pure modes, combined with the use of an average value for σ θθ in several elements ahead of crack, we obtained quite good results.

Numerical application

Several benchmarks were used for testing and validating the numerical implementation of XFEM in DynaCrack. As an illustration of the capabilities of the code, in this paper, we present the numerical simulation of a fracture impact test experimented in our laboratory [START_REF] Nistor | Identification expérimentale et simulation numérique de l'endommagement en dynamique rapide: application aux structures aéronautiques[END_REF]. As shown in Figure 6a, a double-notched specimen is impacted by a cylindrical projectile using a gas-gun device. The shifted angle of the notch induces the crack initiation. The crack propagates during the impact. The propagating direction and final crack length depend on the projectile impact speed. Figure 6b shows a picture of a real specimen impacted at 67 m/s. The material specimen is an aluminium alloy, A-U4G1 T3, with the mass density ρ = 2780 kg/m 3 , Young modulus E = 73.1 GP a and Poisson coefficient ν = 0.33. The cohesive damage crack propagation model parame-ters were identified using a combined Monte-Carlo+Levenberg-Marquardt procedure [START_REF] Nistor | Identification of a Dynamic Viscoplastic Flow Law Using a Combined Levenberg-Marquardt and Monte-Carlo Algorithm[END_REF] as: δ crit = 19.72 • 10 -6 m and σ max = 240 M P a.

The DynaCrack model of this test takes into account its geometrical symmetry and the boundary conditions of the experimental set-up, are shown in Figure 7a. For modeling the impact conditions, the projectile mass was distributed over the impacted face nodes and the impact speed was introduced as initial speed on the same nodes. The mesh accounts 820 elements and a finer meshing was performed for the expected crack zone. Numerical results performed with the XFEM analysis, for experimental corresponding impact speeds which ranges from 40 to 77 m/s, reported a final advancing of the crack from 0.35 to 4.7 mm closely by the experimental ones (0.5 to 4.45 mm). The crack propagation angles ranges from 0˚, for higher impact speed, and 20˚for lower ones. The observed angles of experimental crack propagation paths are about 23˚for higher impact speeds but that are strongly influenced by the experimental setup conditions. The evolution of the shearing field during the analysis for the impact speed of is 67 m/s is illustrated in Figure 7b, c and d. The numerical crack advances over 9 elements, corresponding to a final length of 3.6 mm (the corresponding experimental value is 3.45 mm). The crack path is marked by the red line.

Conclusions

In the present paper, the numerical propagation of a crack is achieved using an explicit XFEM home-made code indented for dynamic analysis. The main challenges of this work concerns the dynamically node enrichment with additional degrees of freedom for the representation of the crack propagation across the elements, the numerical evaluation of the mass and stiffness matrices over the cutted elements, the managing of a variable size for the global degree of freedom vector for post-processing purpose and the implementation of a crack evolution model adapted to our specific crack-tip enrichment. Several algorithms were implemented and tested in order to set-up the DynaCrack module implemented on behalf of the DynELA explicit FEM code for performing XFEM dynamic analysis. A numerical simulation of an impact fracture has been presented and the comparison of numerical and experimental results proves the robustness and the accuracy of the implemented algorithms.

Figure 2 .

 2 Figure 2. Cohesive model description

Figure 3 .

 3 Figure 3. Partition of the cutted elements

Figure 5 .

 5 Figure 5. Cohesive model implementation

Figure 6 .

 6 Figure 6. Design of the impact test used for numerical validation

Figure 7 .

 7 Figure 7. Numerical results for the fracture test model


G Hétreux

R Théry

N Olivier
email: nelly.olivier@ensiacet.fr

J M Le Lann

De la simulation dynamique hybride vers la conduite de procédés batch et semi-continus

Keywords:

PrODHyS est un environnement de simulation dynamique hybride dont la particularité est de fournir des composants orientés objets généraux et réutilisables pour la modélisation d'appareils et d'opérations présents dans les procédés industriels. Les grands principes sur lesquels est fondé PrODHyS sont d'abord exposés en développant plus particulièrement les aspects relatifs au niveau commande. Ces concepts sont ensuite illustrés à travers la simulation d'un exemple académique de procédé batch permettant d'illustrer les potentialités du simulateur. Enfin, la dernière partie aborde brièvement des travaux en cours montrant l'exploitation de ce simulateur pour la conduite de ce type de procédés (ordonnancement et supervision) .

ABSTRACT. PrODHyS is an object-oriented environment dedicated to the hybrid dynamic simulation of industrial processes. The major aspects of the modelling approach proposed by this tool are detailed with a special attention to the command part. Then, these principles are illustrated through the modelling and the simulation of an academic example of batch plant which enables to illustrate the potentialities of the simulator. Finally, the use of this simulator for the scheduling and the monitoring of these processes are described.

MOTS-CLÉS :

Introduction

Aujourd'hui, la simulation dynamique constitue un outil d'étude essentiel pour l'analyse et la conception des procédés industriels. En effet, l'évolution rapide de la puissance des calculateurs et les avancées dans les stratégies de résolution des systèmes d'équations algébro-différentielles, ont permis de valider des approches basées sur l'utilisation de modèles mathématiques rigoureux pour la représentation des opérations unitaires classiques de génie chimique. Dans ce contexte, l'unification des travaux de recherche en modélisation et simulation de procédés menés depuis de nombreuses années au sein du Laboratoire de Génie Chimique a conduit au développement de la plate-forme de simulation dynamique hybride nommée PrODHyS. Si jusqu'à présent, cet environnement a été exploité pour simuler le fonctionnement d'opérations unitaires isolées, les dernières évolutions permettent désormais de traiter des procédés batch complets composés de plusieurs opérations unitaires complexes.

La première partie de cet article permet de présenter et d'illustrer les principes et concepts fondamentaux de PrODHyS par la mise en place de tels systèmes. Pour cela, la section 2 présente brièvement la structure de la bibliothèque et le formalisme hybride mise en oeuvre (Réseaux de Petri Différentiel Objet). La section 3 décrit les principes de modélisation de la partie opérative du procédé, alors que la section 4 s'intéresse plus précisément au niveau commande. Dans la section 5, la modélisation d'un procédé académique permet de mettre en oeuvre l'ensemble de ces éléments. Enfin, dans la section 6, les potentialités de la plateforme sont montrées au travers de travaux de recherche entrepris récemment, visant à exploiter la simulation dynamique hybride pour la conduite de ce type de procédés. Dans ce cadre, une première partie aborde la conception et le développement d'un outil d'ordonnancement fondé sur le couplage de ce simulateur avec un module d'optimisation, alors que la seconde concerne la supervision des procédés et plus précisément les fonctions de surveillance et de diagnostic de fautes en intégrant un modèle de simulation du procédé utilisé comme référence.

L'environnement PrODHyS

Après avoir montré l'intérêt de la simulation dynamique hybride pour l'analyse des procédés batch, les concepts fondamentaux sur lesquels est fondé PrODHyS sont abordés. Dans un second temps, le formalisme hybride mis en oeuvre et la structure générale du modèle de simulation sont brièvement décrits.

Intérêts de la modélisation hybride pour la simulation dynamique des procédés discontinus

Les procédés batch et semi-continus sont souvent les systèmes privilégiés dans les industries du secteur pharmaceutique, cosmétique, de l'agroalimentaire, des polymères ou de la biochimie car ils permettent de fabriquer de petites quantités de produits variés et à haute valeur ajoutée en utilisant les mêmes ressources (équipements, opérateurs, cuves de stockage, matières premières, utilités, etc.). Si par nature, les phénomènes physiques évoluent de manière continue, la modélisation de ces systèmes fait souvent apparaître des aspects événementiels. En effet, la complexité de ces opérations conduit parfois à simplifier ou négliger certains phénomènes physiques. Il en résulte que les modèles dynamiques associés apparaissent alors comme des modèles continus par morceaux (par exemple, changement de phase de la matière) ou même discontinus (par exemple, la commutation d'une vanne Tout ou Rien (TOR)). De plus, ces modèles dynamiques sont souvent à structure variable car le jeu d'équations constituant le modèle de simulation dépend généralement de l'état courant du système. Ceci nécessite donc un noyau capable de gérer automatiquement ces changements de configurations, point qui caractérise spécifiquement les simulateurs dynamiques hybrides [START_REF] Zaytoon | Systèmes dynamiques hybrides[END_REF].

Par ailleurs, ces unités marient généralement des séries d'opérations qui fonctionnent en continu (qu'il s'agisse d'un régime transitoire ou permanent) avec des opérations purement discontinues. De plus, elles sont généralement constituées de nombreux appareils de traitement interconnectés et partagés. Il en résulte que la recette comporte une forte composante événementielle mais dépend aussi de paramètres à valeurs continues (quantités de matière utilisées, conditions opératoires telles que température, pression, concentration, etc.). De ce fait, elle est souvent définie par des événements d'état (seuil de température ou de composition par exemple) et non par des durées ou des dates d'occurrence fixées a priori. C'est pourquoi la simulation de telles opérations et des transformations physico-chimiques associées nécessite la résolution de modèles phénoménologiques hybrides.

Ainsi, les outils classiques de simulation dynamique continue ou de simulation à événements discrets d'ateliers s'avèrent mal adaptés à cette classe de problèmes et la capacité des simulateurs dynamiques hybrides à gérer simultanément et de manière cohérente les comportements continus et discrets devient un atout majeur.

Concepts fondamentaux de l'environnement PrODHyS

Une des bases fondamentales de PrODHyS est la représentation orientée objet des procédés industriels [START_REF] Jourda | Introducing ATOM, the Applied Thermodynamic Object-oriented Model[END_REF][START_REF] Sargousse | Noyau numérique orientée-objet dédié à la simulation des systèmes dynamiques hybrides[END_REF][START_REF] Moyse | Odysseo : plate-forme orientée-objet pour la simulation dynamique des procédés[END_REF], Hétreux et al, 2002). PrODHyS est donc une bibliothèque de classes destinées à être dérivées en exploitant les mécanismes objets (polymorphisme, composition, héritage, généricité). Sa conception repose sur l'identification de composants élémentaires généraux et réutilisables dont l'assemblage permet de construire des modèles de simulation élémentaires, eux-mêmes étant des composants de modèles de simulation plus complexes et spécialisés. La détermination et la conception de ces briques élémentaires constituent un des premiers apports importants de ces travaux.

Le deuxième apport majeur est l'intégration dans PrODHyS d'un noyau de simulation dynamique hybride [START_REF] Hétreux | Object hybrid formalism for modelling and simulation of chemical processes[END_REF], [START_REF] Perret | Intégration des Réseaux de Petri Différentiels à Objets dans une plateforme de simulation dynamique hybride : application aux procédés industriels[END_REF] basé sur le formalisme réseau de Petri différentiel objet (RdPDO).

Architecture logicielle de l'environnement PrODHyS

PrODHyS fournit donc des composants logiciels permettant de traiter des opérations unitaires isolées ou couplées. Leur conception et leur développement s'appuient sur un processus unifié (UP), basé sur le formalisme UML et le langage objet C++.

Elle est organisée en 2 couches fonctionnelles indépendantes (cf. figure 1). La couche simulation regroupe les classes correspondant au noyau de simulation (gestionnaire de simulation, formalisme RdPDO, intégrateurs EDA, solveurs EANL).

Elle fournit donc les objets utiles à la simulation de tout système dynamique hybride, quel que soit le domaine d'application. La couche modélisation rassemble les classes utilisées spécifiquement pour la modélisation des procédés. Actuellement, cette bibliothèque regroupe plus de 1000 classes réparties en 7 modules. Divers outils et IHM [START_REF] Hétreux | Conception d'un atelier de modélisation pour la simulation dynamique hybride des procédés[END_REF] ont été développés (ou sont en cours de développement) selon le contexte d'utilisation de la bibliothèque (figure 1).

Principes et éléments du formalisme RdPDO

La modélisation et la simulation des systèmes hybrides a donné lieu à différents outils et formalismes [START_REF] Zaytoon | Systèmes dynamiques hybrides[END_REF][START_REF] Demongodin | Generalised Batches Petri Net: Hybrid Model for High Speed Systems with Variable Delays[END_REF][START_REF] Champagnat | Modelling and simulation of hybrid system through Pr-Tr PN-DAE model[END_REF][START_REF] Andersson | Object-Oriented Modelling and Simulation of Hybrid Systems[END_REF][START_REF] Barton | The Modelling of Combined Discrete/Continuous Processes[END_REF][START_REF] Nilsson | Dynamic modeling of chemical processes using Omola[END_REF][START_REF] Stephanopoulos | MODEL.LA : a modelling language for process engineering. I-the formal framework[END_REF][START_REF] Piela | ASCEND : an object-oriented computer environment for modelling and analysis[END_REF]. Lors de la conception du noyau de simulation de PrODHyS, il nous a semblé important de le fonder sur un formalisme mixte clairement défini permettant de décrire de manière aussi rigoureuse la dynamique continue que discrète. Si l'utilisation de systèmes d'équations algébro-différentielles paraît incontournable

Pilote de Simulation

Pilote de Simulation

Simulateur général Simulateur général pour la description des aspects continus, l'aspect séquentiel fort des systèmes considérés justifie l'emploi du modèle réseau de Petri. Dans ce cadre, les Réseaux de Petri Différentiel Objet (RdPDO) ont été mis en oeuvre pour spécifier le modèle de simulation associé à chaque composant [START_REF] Hétreux | Object hybrid formalism for modelling and simulation of chemical processes[END_REF]. Comme il n'est pas possible ici de développer de manière exhaustive les fondements théoriques de ce formalisme (voir Perret, 2003, Perret et al, 2004), seuls sont présentés ici les principes élémentaires pour permettre d'appréhender la manière dont les modèles de simulation sont spécifiés via PrODHyS.

Indiquons d'abord que l'approche objet et les réseaux de Petri sont exploités selon une approche dite combinée étendue. Celle-ci consiste à les faire interagir de deux manières : la première est d'introduire les objets dans les réseaux de Petri et la seconde est d'introduire les réseaux de Petri dans les objets. Si la partie statique offre un haut niveau d'abstraction grâce aux concepts objets, la partie dynamique fait collaborer au sein d'une même structure des systèmes d'équations différentielles algébriques (EDA) pour décrire l'évolution continue du système, et des réseaux de Petri utilisés pour spécifier et gérer les séquences légales de commutation entre ces différents systèmes EDA.

Evolution dynamique dans les RdPDO

La figure 2 montre un exemple d'évolution dans les RdPDO. Lorsqu'un jeton marque une place différentielle, le système EDA est instancié par substitution des variables formelles par les attributs du jeton. La résolution induit l'évolution continue de ces attributs ainsi que l'évaluation de la condition associée à la transition (événements d'état ou de temps). Lorsque cette condition est atteinte, l'action est exécutée et la transition est franchie. Plusieurs jetons de même type peuvent marquer la place, déclenchant alors l'intégration du même système EDA à partir de la date de marquage. Par ailleurs, l'instanciation du système EDA d'une place différentielle peut nécessiter une combinaison de jetons de même type et/ou de types différents. Dans ce cas, la cohérence de l'ensemble doit être garantie par la modélisation ou validée a posteriori par la simulation. Notons toutefois qu'il s'agit là d'un problème difficile voire impossible qui reste ouvert. Ainsi, le réseau de Petri peut être vu comme un moniteur d'EDA conduisant à la création dynamique d'un unique modèle de simulation dont la dimension et la structure fluctuent entre deux événements. Outre la résolution des systèmes EDA (basée sur une extension de la méthode de Gear [START_REF] Gear | The simultaneous numerical solution of differential-algebraic equations[END_REF])) et des modèles discrets (joueur spécifique à cette classe de RdP), le noyau effectue d'autres traitements indispensables à ce type de simulateurs tels que le calcul exact des instants de commutation, le recalage des états, la vérification de la cohérence des nouveaux modèles continus générés après commutation, la réinitialisation des variables d'état et de leurs dérivées.

0 64 0 = - 〉 〈 , dt a . y d 0 1 = + 〉 〈 dt a . y d 0 6 = - 〉 〈 a . y 5 ← 〉 〈 a . y {Système EDA 2}
0 64 0 = - 〉 〈 , dt a . y d 0 1 = + 〉 〈 dt a . y d 0 6 = - 〉 〈 a . y 5 ← 〉 〈 a . y {Système EDA 2}

Structure générale du modèle de simulation

Afin de pouvoir réaliser la simulation d'un procédé discontinu, il est nécessaire de modéliser à la fois la partie commande (le superviseur) et la partie opérative (le procédé). Dans PrODHyS (cf. figure 3), le modèle relatif au niveau commande (a priori spécifique à la recette, à la topologie du procédé considéré et au plan de production à réaliser) est complètement dissocié des modèles des appareils qui eux, doivent être réutilisables quelque soit le contexte (notion de composant).

Modélisation de la partie 'opérative'

La modélisation des procédés (appareils et matière qu'ils contiennent) a déjà été présentée dans de précédentes communications [START_REF] Hétreux | Component based approach for simulation of dynamic hybrid systems[END_REF], Perret et al., 2004[START_REF] Perret | Intégration des Réseaux de Petri Différentiels à Objets dans une plateforme de simulation dynamique hybride : application aux procédés industriels[END_REF]. Dans la suite, seuls sont rappelés les principes essentiels de celle-ci. La spécification d'un procédé complet passe d'abord par la spécification de chaque entité élémentaire composant le procédé. Remarquons cependant que ceux-ci sont généralement proposés par les modules Process, Device ou CompositeDevice et évitant ainsi à un utilisateur non modélisateur d'avoir à les spécifier. La spécification de la topologie du procédé complet est alors réalisé par un assemblage plus ou moins complexe de ces objets élémentaires conduisant à la définition de différents niveaux topologiques (cf. figure 3).

Modélisation d'un appareil élémentaire

Tout objet de procédé présent dans PrODHyS est caractérisé par (cf. figure 4) :

• un ensemble de caractéristiques statiques définies au moyen d'une configuration (dimensions, capacité, géométrie de l'appareil, etc.), • un ensemble de variables d'état continues propres à l'entité (température, pression, composition, rétention, débit, etc.)

Modélisation d'un procédé complet

Il s'agit de représenter la structure du procédé selon une vision système, c'est à dire en terme de connexions entre diverses structures composées hiérarchiquement au travers desquelles circulent matière et énergie.

Modélisation topologique selon l'axe connexion

Il s'agit de décrire les connexions entre les différents éléments constitutifs du procédé (figure 3). Une connexion correspond à une visibilité offerte à un élément vers un autre par l'intermédiaire des ports. Elle s'établit par la mise en relation d'un port d'entrée et d'un port de sortie de même type. Toutefois, le modèle reste totalement acausal, la notion d'entrée ou de sortie permettant seulement de fixer le sens conventionnel du flux (positif si sortie entrée). Selon cet axe, un procédé est donc représenté par un graphe orienté.

Modélisation topologique selon l'axe décomposition/composition

L'axe décomposition/composition permet de décrire la structure interne de chaque entité. Fondés sur le design pattern Composite, deux types initiaux d'éléments suffisent à son support : les appareils élémentaires (cf. section 3.1) qui sont non-décomposables et peuvent être spécialisés (fixent la granularité de la modélisation) et les appareils composés qui décrivent récursivement la décomposition. Selon cet axe, un procédé apparaît alors sous la forme d'un arbre dont les feuilles sont les appareils élémentaires et dont la racine est un objet général de type Flowsheet. Par ailleurs, chaque appareil élémentaire issu de la décomposition résulte de la spécialisation de composants encore plus généraux définis dans PrODHyS [START_REF] Hétreux | Bibliothèque orientée-objet pour la simulation dynamique des procédés : architecture et mise en oeuvre[END_REF][START_REF] Moyse | Odysseo : plate-forme orientée-objet pour la simulation dynamique des procédés[END_REF].

Modélisation de la partie 'Commande'

L'évolution d'un procédé discontinu induit une succession de phases transitoires et de régimes continus. De ce fait, un niveau commande est explicitement défini afin de piloter le déroulement de la simulation. Pour cela, cette section décrit la structure du modèle du niveau commande qui doit par ailleurs, assurer la génération quasiautomatique du modèle de simulation à partir de la recette et du plan de production.

F out = 0 F out -F in = 0 F out = d c F out -F in = 0 VANNE V1 ON OFF DETECTEUR D1 V≥T V<T AUTOMATE CUVE ST6 0 F F dt dU out in = + - F in F out F in U V F out EMPTY NOT EMPTY U>0 U=0 Open_V1 Signal_D1 ST6 ∅ = 70 cm h = 110 cm ALIM V1 D1
a -le procédé c -modèle de simu lation (n iveau commande et procédé)

F in F out F in F out U Open_V1 Signal_D1 b -GRAFCET Open_V1 Signal_D1 FILL_ST6 FILL_ST6 ON OFF F out = 0 F out -F in = 0 F out = d c F out -F in = 0 VANNE V1 ON OFF DETECTEUR D1 V≥T V<T AUTOMATE CUVE ST6 0 F F dt dU out in = + - F in F out F in U V F out EMPTY NOT EMPTY U>0 U=0 Open_V1 Signal_D1 ST6 ∅ = 70 cm h = 110 cm ALIM ALIM V1 D1
a -le procédé c -modèle de simu lation (n iveau commande et procédé)

F in F out F in F out U Open_V1 Signal_D1 b -GRAFCET Open_V1 Signal_D1

FILL_ST6

Le marquage d'une place de commande d'un dispositif actif induit en général l'évolution de son RdP, lui-même provoquant éventuellement en cascade l'évolution des dispositifs passifs au travers du réseau constitué par la connexion des différents ports matière ou énergie. Par conséquent, l'évolution des différents modèles est conditionnée par deux types distincts d'événement : Les interactions entre RdP recette et RdP procédé sont illustrées sur la figure 5. Celle-ci représente une séquence opératoire permettant d'alimenter un réacteur jusqu'à ce qu'un volume désiré soit atteint. L'opération de remplissage du réacteur est commandée par la recette qui envoie un jeton sur la place signal de l'appareil vanne. L'alimentation est maintenue ouverte tant que cette place de commande reste marquée. Afin de détecter l'instant où le volume souhaité de réactif est transféré, un détecteur de niveau est utilisé. Le marquage de la place de l'appareil détecteur indique que le volume de réactif a atteint la valeur cible. La transition est alors tirée. L'absence de jeton sur la place commande provoque alors la fermeture de la vanne.

• d'une part,

Représentation hiérarchisée de la recette

La recette regroupe l'ensemble des informations utiles à la fabrication d'un produit. Cependant, lorsque le système de production possède de nombreux équipements ou si l'élaboration d'un produit nécessite de nombreuses phases de transformation, cette recette peut rapidement devenir complexe. C'est pourquoi la norme ISA/SP88 (www.isa.org) propose une modélisation hiérarchisée de la recette et établit une approche standard pour traiter la complexité du contrôle/commande des procédés batch. De manière générale, une recette est constituée de quatre groupes d'information :

• l'en-tête : définit les informations administratives et d'identification,

• la formule : définit les paramètres à appliquer (ingrédients, proportions, conditions opératoires, etc.), • la procédure : définit la suite de tâches ordonnées dans le temps à accomplir pour assurer la fabrication du produit, • les ressources : définit les besoins en équipement pour la fabrication du produit. Comme l'illustre la figure 6, quatre niveaux de recette sont définis, adaptés au niveau de finesse avec lequel le procédé est pris en compte. La recette de contrôle est une recette exécutable qui possède les informations nécessaires à la simulation (volume et composition des lots, équipements utilisés, conditions opératoires, etc.).

Cette structuration de la recette induit aussi une hiérarchisation en cinq niveaux des informations associées à la partie procédurale (figure 7) :

• la procédure définit la séquence d'opérations à exécuter sur une ligne de production pour fabriquer un produit particulier. Une opération correspond à la notion d'opération unitaire induisant une modification des propriétés physicochimiques d'un produit (réaction, séparation, décantation, chauffe, etc.). • chaque opération définit l'ensemble des phases exécutées en séquence ou en parallèle dans un équipement (réacteur, colonne, etc.). La notion de recette principale est représentée comme une spécialisation de la recette de site (MasterRecipe). Elle possède un lien vers la classe ManufacturingCampaign qui regroupe l'ensemble des informations relatives aux tâches impliquées dans la campagne de fabrication. Ces tâches résultent de la prise en compte de la procédure pour chaque OF à réaliser. Cette classe regroupe donc l'ensemble des informations nécessaires à la fonction ordonnancement qui doit effectuer un découpage de chaque OF en batch (détermination de leur nombre et de leur taille) compatibles avec la capacité de chaque ressource, et calculer les dates de lancement de chacun de ces lots compte tenu de la procédure et de la disponibilité des ressources.

Enfin, la recette de contrôle regroupe les informations utiles à la réalisation des lots et au pilotage des équipements. Il s'agit notamment des conditions opératoires (température, pression, composition, etc.) utilisées pour définir les fonctions de transition sur événements d'état au sein des RdPDO du niveau commande. Ici encore, la décomposition présentée précédemment a été adoptée pour la modélisation de la procédure sous PrODHyS. Dans ce cadre, la notion de macroplace a été introduite. Elle permet de remplacer une séquence de places/transitions relatives à une opération ou à une phase par une unique macroplace (figure 10). Cette séquence est alors délimitée par 2 places particulières E et S entre lesquelles les différents types de places définis dans les RdPDO peuvent apparaître (y compris d'autres macro-places). On retrouve ainsi la classe ProcedureRecipe définie comme une séquence d'opérations (classe OperationRecipe), elles-mêmes définies comme un ensemble de phases (classe PhaseRecipe). Une phase est alors définie par un RdPDO sans macro-place, où chaque place correspond à la notion de pas. La notion d'instruction est, dans notre cas, souvent confondue avec la notion de pas.

Application à un procédé didactique

Les différents éléments présentés dans les sections précédentes sont illustrés au travers de la modélisation et la simulation d'un procédé académique, typique des installations batch trouvées en milieu industriel.

Caractéristiques de la recette associée au procédé considéré

Formulation

Le but de l'installation est de fabriquer et de conditionner un produit P dont la pureté doit être au moins égale à 98 % molaire. La réaction considérée est une réaction équilibrée endothermique impliquant deux réactifs R1 et R2 et le produit P et qui suit une loi cinétique d'ordre 2 :

(R) : R1 + R2 ↔ P P RT Ea 1 2 R 1 R RT Ea 1 R x . e . k x . x . e . k r 1 1 - - °- - °- =
Comme le montre la figure 11, la fraction molaire de produit P ne peut excéder 0,83 et un temps de séjour infini serait nécessaire pour atteindre cette composition. Afin de maximiser le taux de conversion de la réaction sans pour autant pénaliser le temps de cycle du procédé, la réaction est stoppée dés que la composition en produit P atteint la valeur de 0,8. Par ailleurs, compte tenu de la loi de vitesse considérée, la réaction (R) est d'autant plus rapide que la température T est élevée. La température retenue pour la réaction doit garantir une réaction rapide tout en maintenant les constituants dans l'état liquide. Une température de 383 K permet de satisfaire ces 2 contraintes. La spécification de pureté imposée pour le produit P étant égale à 98%, une étape de distillation est donc nécessaire. o laisser réagir jusqu'à ce que la composition en produit P atteigne la valeur de 0,8. o vidanger tout ou partie du mélange réactionnel. Lorsque cette opération est effectuée dans le réacteur BR2, les réactifs peuvent aussi provenir du bac de stockage ST3 qui contient un mélange de produit R1 et R2 proche des proportions stoechiométriques issu du recyclage provenant de l'étape de séparation.

Stockage du mélange après réaction dans ST5 ou ST6 selon leur disponibilité

Purification du produit avec le séparateur SEP

3.1 Chargement du rebouilleur avec une partie du contenu de la cuve ST5.

Remplissage et démarrage de la colonne 3.3. Mise en régime de la colonne

o transférer en continu le mélange du réservoir ST5 sur le plateau d'alimentation de la colonne à un débit constant (pompe P4). Ceci provoque le remplissage du rebouilleur o réaliser l'opération de distillation à reflux infini pendant 1h.

Marche à reflux fini

Ouvrir le diviseur de courant VPA1 afin de réaliser l'opération de distillation à reflux fini (r=1.2). Le produit P quasiment pur est récupéré dans le réservoir ST4. Le mélange de réactifs résiduels R1 et R2 est recyclé dans le réservoir ST3.

Stockage des réactifs n'ayant pas réagi dans le réservoir ST3

Ceux-ci sont recyclés dans le procédé et peuvent servir de matières premières à la réaction mise en oeuvre dans le réacteur BR2. Lorsque cette opération est effectuée dans le réacteur BR2, les réactifs peuvent aussi provenir du bac de stockage ST3 qui contient un mélange de produit R1 et R2 proche des proportions stoechiométriques issu du recyclage provenant de l'étape de séparation.

Conditionnement du produit P

Stockage du mélange après réaction dans ST5 ou ST6 selon leur disponibilité

Purification du produit avec le séparateur SEP

Chargement du rebouilleur avec une partie du contenu de la cuve ST5.

Remplissage et démarrage de la colonne 3.3. Mise en régime de la colonne

o transférer en continu le mélange du réservoir ST5 sur le plateau d'alimentation de la colonne à un débit constant (pompe P4). Ceci provoque le remplissage du rebouilleur o réaliser l'opération de distillation à reflux infini pendant 1h.

Marche à reflux fini

Ouvrir le diviseur de courant VPA1 afin de réaliser l'opération de distillation à reflux fini (r=1.2). Le produit P quasiment pur est récupéré dans le réservoir ST4. Le mélange de réactifs résiduels R1 et R2 est recyclé dans le réservoir ST3.

Stockage des réactifs n'ayant pas réagi dans le réservoir ST3

Ceux-ci sont recyclés dans le procédé et peuvent servir de matières premières à la réaction mise en oeuvre dans le réacteur BR2.

Conditionnement du produit P

Le produit P, stocké dans le réservoir ST4 est ensuite conditionné dans des bidons de 10 litres

Fabrication du produit P selon la réaction (R) dans le réacteur BR1 ou BR2,

Cette opération comporte les phases suivantes : o charger le réacteur avec n/2 moles de produit R1, o préchauffer à 383 K, o charger le réacteur de n/2 moles de produit R2 tout en régulant la température à 383 K. o laisser réagir jusqu'à ce que la composition en produit P atteigne la valeur de 0,8. o vidanger tout ou partie du mélange réactionnel. Lorsque cette opération est effectuée dans le réacteur BR2, les réactifs peuvent aussi provenir du bac de stockage ST3 qui contient un mélange de produit R1 et R2 proche des proportions stoechiométriques issu du recyclage provenant de l'étape de séparation.

Stockage du mélange après réaction dans ST5 ou ST6 selon leur disponibilité

Purification du produit avec le séparateur SEP

Chargement du rebouilleur avec une partie du contenu de la cuve ST5.

Remplissage et démarrage de la colonne 3.3. Mise en régime de la colonne

o transférer en continu le mélange du réservoir ST5 sur le plateau d'alimentation de la colonne à un débit constant (pompe P4). Ceci provoque le remplissage du rebouilleur o réaliser l'opération de distillation à reflux infini pendant 1h.

Marche à reflux fini

Ouvrir le diviseur de courant VPA1 afin de réaliser l'opération de distillation à reflux fini (r=1.2). Le produit P quasiment pur est récupéré dans le réservoir ST4. Le mélange de réactifs résiduels R1 et R2 est recyclé dans le réservoir ST3.

Stockage des réactifs n'ayant pas réagi dans le réservoir ST3

Ceux-ci sont recyclés dans le procédé et peuvent servir de matières premières à la réaction mise en oeuvre dans le réacteur BR2.

Conditionnement du produit P

Modélisation de la partie opérative

A titre d'illustration, la figure 14 montre la structure topologique interne de l'objet colonne selon l'axe composition/décomposition. Une description plus complète de la modélisation de cet objet peut être trouvée dans [START_REF] Théry | Environnement orienté objet pour la simulation dynamique d'un procédé de distillation réactive[END_REF]. Par ailleurs, la modélisation des autres équipements composant ce procédé est donnée dans [START_REF] Perret | Intégration des Réseaux de Petri Différentiels à Objets dans une plateforme de simulation dynamique hybride : application aux procédés industriels[END_REF].

Diviseur de courant

Diviseur de courant

Diviseur de courant

Diviseur de courant

Vanne

Vanne Vanne Vanne

Modélisation de la partie commande

Il n'est pas possible de donner ici le modèle complet du niveau commande associé à ce procédé. Aussi, seuls les RdPDO découlant de la décomposition de la procédure selon les différents niveaux qui correspondent à l'opération de séparation sont décrits dans la suite.

Le niveau procédure

La figure 15 montre le modèle associé à la procédure (spécifique à un mode de production fixé). Celle-ci est constituée des 4 opérations unitaires, chacune étant représentée par une macro-place : synthèse dans BR1, synthèse dans BR2, séparation, conditionnement et gestion des cuves de stockage ST5 et ST6. Dans ce cadre, le lien prévu vers un objet de type MasterRecipe trouve ici sa justification car l'établissement de la procédure requiert la connaissance de l'ensemble des informations utiles pour la fabrication du produit (gammes de fabrication, équipements nécessaires).

Le niveau opération

Ce niveau de description de la procédure rejoint la notion d'opération unitaire couramment rencontrée en génie des procédés. La figure 16 représente l'opération de séparation, associée à la macro-place M30 dans le RdP de la procédure. L'expansion de M30 se compose d'une phase composée de distillation représentée par une macro-place notée M300 et de phases élémentaires (non décomposables) représentées par une place (place nommée FILL_REB, par exemple).

Le niveau phase

Une phase est une succession de pas permettant de réaliser une fonction élémentaire. La figure 17 représente l'expansion la macro-place M300 associée à la phase de distillation (cf. figure 16). Le RdPDO représentant le fonctionnement intrinsèque de la colonne est aussi représenté. Comme nous pouvons le constater les évolutions respectives de ces deux RdP sont intimement liées : • l'évolution du RdP associé à la phase M300 active des places de commande permettant l'évolution du RdP intrinsèque de la colonne.

• ce dernier permet lui-même de faire évoluer les RdP représentant l'évolution des objets composant la colonne (ici, le condenseur et le rebouilleur). • inversement, l'évolution du RdP intrinsèque de la colonne active des places d'information qui permettent l'évolution du RdP de la phase de distillation.

Précisons par ailleurs, qu'au niveau de la phase, les RdP sont définis de manière générique afin de permettre leur réutilisation dans un autre contexte. Les modèles continus associés à chaque objet constituant ce procédé prennent en compte, entre autres, les bilans matière globaux et partiels, le bilan d'énergie, les équilibres liquide/vapeur, les vitesses de réaction, ainsi que les phénomènes hydrauliques. En effet, hormis les conduites munies d'une pompe, les transferts entre cuves sont effectués par gravité. A titre d'illustration, la figure 18 donne une partie de l'évolution de composition au sein du réacteur BR1 durant l'opération de réaction. La capacité du simulateur à gérer les évènements d'état est un atout majeur car cela permet de :

• prendre en compte des événements dont la date d'occurrence n'est pas connue a priori, • représenter avec plus de précision le comportement du procédé puisqu'il s'appuie sur la simulation des phénomènes physiques via des modèles prenant en compte les conditions opératoires (température, pression, etc.)

En effet, l'hypothèse d'une durée constante et connue des tâches constitue une restriction sévère face à la sensibilité de certaines opérations unitaires à ces paramètres. La colonne à distiller discontinue est un autre exemple où la durée opératoire est liée à plusieurs facteurs (qualité de la charge initiale, politique de chauffe au bouilleur, reflux en tête de colonne, débits de soutirage latéraux, pertes thermiques, etc). De plus, la durée d'une tâche peut aussi dépendre de l'état du système à un instant donné. Par exemple, la durée de chauffe d'un produit dépend de sa température initiale, elle-même pouvant dépendre de la durée d'attente du produit dans la cuve de stockage amont si des pertes thermiques existent. De même, lorsque les transferts entre cuves sont opérés par gravité, le temps de transfert dépend du débit qui devient une fonction du niveau de produit dans la cuve source à cause des phénomènes de pression hydraulique. L'estimation de cette durée est rendue encore plus délicate si des alimentations et des vidanges sont autorisées simultanément. Ne pas inclure ces paramètres dans la représentation du système impose une approximation de la durée opératoire. Afin d'évaluer l'impact de cet aspect, une série de simulations conduites successivement sur événements de temps et sur événements d'état a été réalisée pour un procédé comportant des phases de transfert par gravité entre cuves de même géométrie. Les durées de transfert ont été évaluées en supposant successivement que le débit est minimum (fond de cuve), maximum (cuve pleine) ou moyen (milieu de cuve). Les résultats obtenus pour les différents scenarii envisagés sont reportés dans le tableau 1. Les écarts sur les temps de cycle entre les 2 types de simulation sont compris entre 3,5% (débit minimum) et 14,6% (débit maximum). En considérant le débit moyen, l'écart est de 5,6%. On constate donc que les écarts ne sont pas toujours négligeables et peuvent nuire à la qualité des résultats obtenus, notamment si l'on cherche à optimiser le fonctionnement de l'unité. De façon générale, la politique de transfert des lots entre appareils est un facteur important de productivité d'une unité batch. En effet, la possibilité de pouvoir accumuler des produits intermédiaires dans des cuves tampon de capacité finie (telles que ST1, ST2, ST3, ST4, ST5 et ST6) offre un maximum de souplesse dans la gestion des lots en découplant les phases de vidange et de remplissage de deux opérations de transformation successives. Une cuve tampon permet aussi de modifier le volume du lot entre les opérations amont et aval, avec un choix de volume appartenant à une gamme continue. Or, la taille et le nombre de lots nécessaire pour satisfaire un objectif de production est souvent une donnée d'entrée car ces paramètres sont difficiles à établir seulement par simulation.

Choix dans l'affectation des ressources

Dans un atelier où plusieurs ressources peuvent réaliser la même tâche, le choix dans l'affectation d'une ressource peut influencer sensiblement les performances de l'atelier. Or, choisir d'utiliser le réacteur BR1 ou BR2 à un instant donné, par exemple, peut dépendre de la disponibilité des matières premières dans les cuves de stockage amont (ST1, ST2, ST3), mais aussi de la capacité résiduelle des cuves de stockage aval (ST5 et ST6). Là encore, l'affectation des tâches aux ressources est souvent une donnée d'entrée car la simulation dispose de peu de moyen pour faire ce choix judicieusement.

Transition entre mode de production

La présence d'un réservoir tampon représente la configuration la plus classique pour assurer la transition entre deux zones de procédé ayant un mode de production différent. Cependant, la gestion de ces transitions est particulièrement délicate en cours de simulation.

En effet, lorsque la matière circule de la zone discontinue vers la zone continue du procédé, l'activité de la zone aval est conditionnée à l'existence d'une quantité suffisante de matière dans le réservoir. Si le débit de soutirage de la zone continue provoque une baisse de niveau permanente du réservoir qui n'est pas compensée par l'arrivée des lots de la zone amont, une rupture de production de la zone continue est inéluctable au bout d'un certain temps. Dans notre exemple, les cuves ST5 et ST6 jouent précisément ce rôle de régulateur en assurant un flux continu de matière entre la zone réactionnelle qui fonctionne en mode batch et la colonne de distillation qui fonctionne en continu, une fois démarrée. La gestion de telles situations est particulièrement importante dans la mesure où la remise en régime d'une installation continue (telle qu'une colonne de distillation) est souvent longue et consommatrice d'énergie. A l'inverse, si le soutirage continu est trop faible, il y aura dépassement de la capacité de stockage et l'activité de l'atelier discontinu devra être freinée. Cette situation qui semble plus facile à gérer peut toutefois être délicate lorsque les lots ne peuvent patienter dans les appareils que durant un temps limité. Dans ces configurations, la détermination puis le respect des dates de lancement d'opération dans la zone discontinue devient une contrainte prépondérante.

Pour un sens de circulation allant du continu vers le discontinu, la gestion du réservoir tampon obéit aux mêmes principes. En effet, la rupture de production du procédé continu est déclenchée par l'atteinte d'une valeur seuil maximale de la capacité du réservoir. Inversement, l'absence de quantité de matière suffisante dans le réservoir pour lancer la fabrication d'un lot retarde son démarrage et diminue la productivité de la chaîne discontinue. Ceci correspond typiquement aux situations qui doivent être gérées au niveau des cuves ST3 et ST4.

Exploitation de la simulation dynamique pour la conduite des procédés

Cette dernière section présente des travaux en cours exploitant le simulateur PrODHyS. La première application concerne la conception et le développement d'un outil d'ordonnancement basé sur le couplage de ce simulateur avec un module d'optimisation de type MILP. La deuxième application intègre le simulateur au sein d'un système de supervision afin d'assurer les fonctions de surveillance et de diagnostic.

6.1. Vers un ordonnancement réactif des procédés 6.1.1. Motivations du couplage simulation/optimisation pour l'ordonnancement Gérer les problèmes soulevés dans la section 5.4 en cours de simulation apparaît donc difficile. L'idée est donc de coupler la simulation avec un module d'optimisation s'appuyant sur un modèle mathématique plus apte à résoudre ces situations. Le rôle de ce dernier est donc de fournir les éléments de décisions qui manquent à la simulation pour arriver à son terme et non d'explorer l'ensemble des solutions en vue d'établir un ordonnancement « optimal ».

Cependant, dans un contexte plus général, l'exploitation rationnelle de la flexibilité intrinsèque des procédés batch confère à la fonction ordonnancement un rôle crucial. De nombreux travaux dans la littérature traitent ce problème reconnu complexe et différentes approches ont été proposées. Parmi celles-ci, les méthodes d'optimisation basées sur des formulations MILP (Mixed Integer Linear Programming) ou MINLP (Mixed Integer Non Linear Programming) sont souvent mises en oeuvre. Cependant, comme l'indique [START_REF] Mendez | State of the art review of optimisation methods for short term scheduling of batch processes[END_REF], il subsiste un écart entre théorie et pratique, soit à cause de la taille des systèmes considérés, soit à cause des hypothèses simplificatrices qu'il faut éventuellement introduire pour rendre le problème solvable (l'une et l'autre étant d'ailleurs souvent liées). D'un autre coté, de nombreux facteurs induisent que la recherche d'un optimum mathématique sur des modèles trop simplifiés parait souvent inutile dans la pratique. En effet, l'optimalité d'une telle solution est souvent perdue à cause de la nature dynamique de ces systèmes industriels (aléas, évolution de la demande, pannes, etc.). Par ailleurs, la mise en oeuvre d'une solution est rarement directe dés lors qu'elle est confrontée aux caractéristiques du procédé réel (durées opératoires approximées, hypothèses simplificatrices, etc.).

Lorsque la garantie d'optimalité n'est plus le critère central, certaines approches proposées consistent à établir une solution initiale a priori « sous-optimale » (par rapport à certains critères) que l'on tente ensuite d'améliorer progressivement et ponctuellement par différentes techniques (Méndez et Cerda, 2003), [START_REF] Roslöf | An MILP-based reordering algorithm for complex industrial scheduling and rescheduling[END_REF]. La solution initiale est obtenue dans un temps fixé ou pour un effort de calcul donné, soit par une méthode exacte, soit par une méta-heuristique. D'autres travaux sont fondés sur le couplage d'un algorithme génétique avec un simulateur à événements discrets [START_REF] Baudet | A discrete-event simulation approach for scheduling of batch processes[END_REF]. ). Cette procédure itérative est effectuée jusqu'à ce que l'opérateur estime la solution satisfaisante. Par ailleurs, la reproductibilité du séquencement des lots sur les équipements entre deux phases d'optimisation peut constituer un indicateur de convergence de la méthode. Le principe fondamental de l'approche est donc de supposer que la solution « approchée » (en terme de comportement) établie par le modèle d'optimisation avec un effort de calcul réduit est compensée par la modélisation plus fine du procédé réalisé au niveau simulation, notamment en prenant en compte des modes de fonctionnement particuliers (chevauchement des transferts, par exemple). Par ailleurs, cette approche devrait globalement rendre plus robuste les plans de production établis et en facilite l'analyse physico-chimique des phénomènes qui s'y déroulent.

Gestion des alarmes non non

Contraintes violées

Gestion des alarmes non non

Contraintes violées

Gestion des alarmes non non

Contraintes violées

Responsable

Compte tenu de l'ampleur du travail nécessaire pour automatiser cette procédure, ce système est encore en cours de développement. Toutefois, la formulation MILP du modèle d'optimisation est décrite dans la section suivante.

Brève présentation du module d'optimisation ProSched

Diverses formulations MILP pour le problème d'ordonnancement des procédés batch peuvent être trouvées dans la littérature [START_REF] Shah | Integrated batch plant design: a polymer plant case study[END_REF], [START_REF] Floudas | Continuous-time versus discrete-time approaches for scheduling of chemical processes : a review[END_REF], [START_REF] Burkard | Review, extensions and computational comparison of MILP formulations for scheduling of batch processes[END_REF], [START_REF] Shaik | Continuous-Time Models for Short-Term Scheduling in Multipurpose Batch Plants : a Comparative Study[END_REF], (Méndez et al, 2006). Compte tenu de la procédure à mettre en oeuvre, les formulations MILP en temps continu apparaissent les mieux adaptées à nos besoins. Une comparaison détaillée de ce type de modèle peut être consultée dans [START_REF] Shaik | Continuous-Time Models for Short-Term Scheduling in Multipurpose Batch Plants : a Comparative Study[END_REF]. Plus précisément, une formulation du type Unit-Specific Event-Based proposée dans [START_REF] Ierapetritou | Effective continuous-time formulation for short-term scheduling: 1. Multipurpouse batch processes[END_REF] L'ordonnancement obtenu par la résolution de ce modèle permet d'obtenir entre autres, les dates de lancement des lots au niveau de la simulation qui évitent une rupture d'alimentation de la colonne en régime continu. Notre premier objectif (mené à son terme la simulation en respectant toutes les contraintes de production) est donc atteint [START_REF] Hétreux | Approche combinée Optimisation/Simulation pour l'ordonnancement des procédés discontinus[END_REF]. En revanche, la phase de rebouclage entre simulation et module d'optimisation (démarche d'amélioration) n'a pas encore été implementée actuellement au sein de notre système et constitue le prochain objectif de ces travaux. Le module de surveillance doit être capable de traiter aussi bien les défauts du système physique (fuites, pertes thermiques, etc.) que ceux des appareils de contrôle/commande (actionneurs, capteurs, etc.). Comme définie dans (De [START_REF] Kleer | Diagnosing multiple faults[END_REF], notre approche est orientée « bon fonctionnement ». Elle repose donc sur la comparaison du comportement attendu, obtenu à l'aide de la simulation du modèle de référence (valeurs des variables d'état) et le comportement réel observé (mesures provenant du procédé). Par ailleurs, une détection par seuil fixe est utilisée. Cette méthode consiste à détecter l'instant pour lequel les grandeurs surveillées sortent d'une zone délimitée par des seuils autour d'une valeur objectif (valeur de la grandeur au sein du modèle de référence). Bien que des méthodes telles que la détection par seuil variable, soient plus rigoureuses et permettent la minimisation du nombre de détections erronées, [START_REF] Duviella | Conduite réactive des systèmes dynamiques étendus à retards variables[END_REF] indique que la détection par seuil fixe reste plus robuste aux perturbations. Par ailleurs, elle permet d'obtenir une décision assez rapidement et est facilement implémentable et paramétrable. Pour assurer un fonctionnement cohérent de cette fonction, il faut d'abord que les mesures soient filtrées afin d'éliminer le bruit. Toutefois, il peut subsister une incertitude spatiale et l'écart ε entre mesures et variables du modèle au-delà duquel on considère qu'un défaut est avéré, reste un point délicat à évaluer. Généralement, cette valeur ε est un compromis obtenu à partir d'une série de simulations dans lesquelles sa valeur est ajustée. D'autre part, la simulation du modèle de référence doit évoluer en temps réel et de manière synchrone avec le procédé surveillé. En général, cette caractéristique ne peut pas être garantie à tout instant. En effet, une incertitude temporelle existe du fait des paramètres intrinsèques des modèles de simulation décrivant la dynamique du système. Ainsi, certains événements temporels ou d'état peuvent apparaître au niveau du modèle de référence simulé soit en avance, soit en retard par rapport aux observations faites sur le procédé surveillé. Ici encore, se pose le problème de déterminer à partir de quel écart ε peut-on considérer qu'une défaillance est potentiellement détectée. Dans tous les cas, une procédure de recalage du modèle de référence sur le procédé réel doit être réalisée pour pouvoir valider un test de détection.

Lorsque l'occurrence d'un défaut est avérée, la phase de diagnostic est lancée. Ce module s'appuie sur un outil intégrant des connaissances qualitatives du procédé dans lequel sont répertoriées des règles du type « SI condition ALORS cause » (méthodes de classification [START_REF] Dubuisson | Diagnostic, intelligence artificielle et reconnaissance de formes[END_REF]). Toutefois, pour affiner la discrimination et l'interprétation des défaillances constatées, il parait intéressant d'exploiter les informations fournies par la simulation du modèle de référence, notamment la valeur des variables d'état ainsi que leurs dérivées. En effet, si les mesures provenant du procédé sont pertinentes pour la commande, elles sont souvent insuffisantes pour le diagnostic. Des informations complémentaires concernant ces travaux peuvent être trouvés dans (Olivier et al, 2006a), (Olivier et al, 2006b).

Figure 1 .

 1 Figure 1. Environnement PrODHyS et outils associés

Figure 2 .

 2 Figure 2. Principe d'évolution des RdPDO

Figure 3 .

 3 Figure 3. Modélisation topologique d'un procédé de distillation (niveau 1)

Figure 4 .

 4 Figure 4. Exemple simplifié de modélisation d'un objet de procédé élémentaire

4. 1 .

 1 Connexion entre la partie commande et la partie opérative Le modèle associé au niveau commande est un RdPDO maître dont l'évolution provoque l'évolution des RdPDO des entités du niveau procédé. Ce RdPDO, appelé RdP Recette, correspond en quelque sorte au programme GRAFCET qui serait exécuté par un automate programmable industriel (API) pour réaliser un cycle de travail (cf. figure 5). Les signaux échangés entre la partie procédé et la partie commande correspondent soit à l'émission d'une commande (signal Open_V1), soit à la réception d'une information provenant d'un capteur (signal Signal_D1). Un signal est modélisé par une place respectivement dite de commande ou d'information et son état est associé au marquage de cette place. Ces places sont le seul lien entre les appareils commandés et la recette. Considérés comme des boites noires, on distingue alors deux types de dispositif : • les dispositifs actifs : appareils dont le RdP possède une ou plusieurs places de commande et/ou d'information tels que les vannes, les pompes, les alimentations d'énergie, les capteurs, etc (cf. VANNE V1, DETECTEUR D1), • les dispositifs passifs : entités dont le RdP ne possède aucun lien direct avec le RdP recette tels que les cuves, les réacteurs ou la matière (cf. CUVE ST6).

Figure 5 .

 5 Figure 5. Interaction entre le niveau commande et le niveau procédé

Figure 6 .

 6 Figure 6. Hiérarchisation de la recette selon la norme ISA-SP88

Figure 7 .

 7 Figure 7. Décomposition de la procédure selon la norme ISA-SP88

Figure 8 .

 8 Figure 8. Extrait du diagramme de classes associé à la notion de recette

Figure 9 .

 9 Figure 9. Extrait du diagramme de classes associé à la notion de procédure

Figure 10 .

 10 Figure 10. Macro-place sous PrODHyS

Figure 11 .

 11 Figure 11. Vitesse de la réaction en fonction de la température

Figure 12 .Figure 13 .

 1213 Figure 12. Procédure de fabrication du produit P

Figure 14 .

 14 Figure 14. Structure de l'objet colonne (niveaux topologiques 2, 3 et 4)

Figure 15 .

 15 Figure 15. Réseau de Petri décrivant la procédure

Figure 16 .

 16 Figure 16. Expansion de la macro-place associée à l'opération de séparation

Figure 17 .

 17 Figure 17. Expansion de la macro-place M300 représentant la phase de distillation et connexion avec les places de commande et d'information de l'objet colonne

Figure 18 .

 18 Figure 18. Evolution de la composition dans le réacteur BR1

  Comparaison entre simulation sur événements temporel et événements d'état 5.4.2. Détermination de la taille des lots et politique de transfert des lots

Tableau 2 .Inconvénients•

 2 Optimisation/Simulation pour l'ordonnancement Ordonnancement basé sur l'optimisation Ordonnancement basé sur la simulation Avantages • Exploration exhaustive des solutions candidates • Prise en compte globale des contraintes • Méthodes efficaces de résolution • Modélisation souvent plus fidèle du procédé réel • Durées opératoires évaluées par des modèles phénoménologiques Modélisation souvent basée sur des hypothèses simplificatrice ne permettant plus d'exploiter toute la flexibilité de l'atelier • Durées opératoires fixées et souvent surestimées • Evaluation d'une solution candidate simulation du procédé pour un séquencement et une taille de lot donnés Exploration difficile et coûteuse en temps CPU • Difficultés à prendre en compte les contraintes de date (nowait, calendrier au conditionnement, politique de nettoyage) car vision partielle de l'horizon temporel Pour notre part, l'outil d'ordonnancement proposé s'appuie sur ce principe d'amélioration à travers le couplage d'un module d'optimisation basé sur une modélisation de type MILP avec notre simulateur dynamique hybride. Le tableau 2 récapitule les principaux points forts et points faibles de ces deux outils que leur association peut compenser. 6.1.2 Processus de décision et structure de l'outil d'ordonnancement proposé L'établissement d'un plan est basé sur un processus itératif entre le module d'optimisation ProSched et le simulateur PrODHyS dans lequel l'opérateur est intégré (cf. figure 19).

Figure 19 .

 19 Figure 19. Structure de l'outil d'ordonnancement

Figure 20 .

 20 Figure 20. Variables de décision dans la formulation Unit-Specific Event-Based

Figure 22 .

 22 Figure 21. STN (State-Task Network) du procédé considéré

MODELISATION MODELISATION SIMULATION SIMULATION Composite device

	Appareils obtenus
	par spécialisation
	et composition
	de classes du
	module Device

Reaction

Modélisation objet des réactions chimiques

Device Appareils élémentaires concrets construits avec les classes des modules Process, ATOM et Reaction ATOM Base de données thermodynamiques et représentation objet de la matière DISCo Noyau numérique et représentation objet des modèles continus

Hybrid Formalisme RdPDO et noyau de simulation hybride

	MODELISATION MODELISATION
	Process
	Description
	générale des
	procédés
	Composite
	device
	Appareils obtenus
	par spécialisation
	et composition
	de classes du
	module Device
	SIMULATION SIMULATION

Reaction

Modélisation objet des réactions chimiques

Device Appareils élémentaires concrets construits avec les classes des modules Process, ATOM et Reaction ATOM Base de données thermodynamiques et représentation objet de la matière DISCo Noyau numérique et représentation objet des modèles continus

Hybrid Formalisme RdPDO et noyau de simulation hybride

		Process
	Générateurs de code	Description générale des	PrODHyS WinTester
		procédés
	Atelier de Modélisation
			Simulateurs dédiés

IHM Windows pour la saisie du procédé IHM Windows pour la saisie du procédé Pilote de Simulation Pilote de Simulation Simulateur général Simulateur général MODELISATION MODELISATION SIMULATION SIMULATION Composite device

	Appareils obtenus
	par spécialisation
	et composition
	de classes du
	module Device

Reaction

Modélisation objet des réactions chimiques

Device Appareils élémentaires concrets construits avec les classes des modules Process, ATOM et Reaction ATOM Base de données thermodynamiques et représentation objet de la matière DISCo Noyau numérique et représentation objet des modèles continus

Hybrid Formalisme RdPDO et noyau de simulation hybride

	MODELISATION MODELISATION
	Process
	Description
	générale des
	procédés
	Composite
	device
	Appareils obtenus
	par spécialisation
	et composition
	de classes du
	module Device
	SIMULATION SIMULATION

Reaction

Modélisation objet des réactions chimiques

Device Appareils élémentaires concrets construits avec les classes des modules Process, ATOM et Reaction ATOM Base de données thermodynamiques et représentation objet de la matière DISCo Noyau numérique et représentation objet des modèles continus

Hybrid Formalisme RdPDO et noyau de simulation hybride

		Process
	Générateurs de code	Description générale des	PrODHyS WinTester
		procédés
	Atelier de Modélisation
			Simulateurs dédiés

IHM Windows pour la saisie du procédé IHM Windows pour la saisie du procédé

  des événements dits externes qui provoquent des commutations commandées. Il s'agit de tous les signaux échangés entre le niveau commande et le niveau procédé (soit les commandes émises depuis le RdP recette afin de piloter les dispositifs actifs, soit à l'occurrence d'un événement d'état (détection d'un seuil) ou temporel. Spécifiés par l'utilisateur, ces événements apparaissent explicitement sur le RdP recette.

• d'autre part, des événements dits intrinsèques dont l'occurrence dépend uniquement de l'évolution spontanée du procédé. Ces commutations autonomes correspondent au changement d'état de la matière (passage de l'état liquide à l'état liquide/vapeur lorsque la température d'ébullition est atteinte, par exemple) ou d'un dispositif passif. Ceux-ci n'apparaissent donc pas explicitement sur le RdP de la recette (l'utilisateur n'a donc pas à les spécifier) et sont traités exclusivement au sein du modèle de l'entité concernée.

Formule En-tête Procédure Générique Générique Spécifique Spécifique

  Recette adaptée au site, elle précise les caractéristiques des équipements et des ingrédients requis. Elle est utilisée au niveau de la planification.

		Besoins en équipements Besoins en équipements	est une partie de est une partie de	est une partie de est une partie de	En-tête En-tête
		Formule Formule	est une partie de est une partie de	est une partie de est une partie de	Procédure Procédure
		Générique Générique
	Spécifie la méthode générique d'obtention du produit Spécifie la méthode générique d'obtention du produit Spécifie la méthode générique d'obtention du produit Spécifie la méthode générique d'obtention du produit
	(matériaux, proportions, protocole). A ce niveau, elle est (matériaux, proportions, protocole). A ce niveau, elle est (matériaux, proportions, protocole). A ce niveau, elle est (matériaux, proportions, protocole). A ce niveau, elle est
	indépendante d'un site particulier (installation), ou d'une technologie particulière (équipements) et la formule du produit est définie par des proportions. indépendante d'un site particulier (installation), ou d'une technologie particulière (équipements) et la formule du produit est définie par des proportions. indépendante d'un site particulier (installation), ou d'une indépendante d'un site particulier (installation), ou d'une technologie particulière (équipements) et la formule du technologie particulière (équipements) et la formule du produit est définie par des proportions. produit est définie par des proportions.	Est utilisée pour mettre à jour	Est transformée en Est transformée en
					Informations spéci fiques Informations spéci fiques Informations spéci fiques
	Recette adaptée au site, elle précise les caractéristiques Recette adaptée au site, elle précise les caractéristiques Recette adaptée au site, elle précise les caractéristiques				au si te au si te au si te
	des équipements et des ingrédients requis. Elle est des équipements et des ingrédients requis. Elle est des équipements et des ingrédients requis. Elle est
	utilisée au niveau de la planification. utilisée au niveau de la planification. utilisée au niveau de la planification.	Est utilisée pour mettre à jour mettre à jour Est utilisée pour	Est transformée en Est transformée en
					Informations spéci fiques Informations spéci fiques Informations spéci fiques
					au jeu d'équipements au jeu d'équipements au jeu d'équipements
		Est utilisée pour mettre à jour Est utilisée pour mettre à jour	Est transformée en Est transformée en
					Informations relatives aux Informations relatives aux
					lots, opérateurs, aux lots, opérateurs, aux
		Spécifique Spécifique			ressources … ressources …

Recette de site adaptée à une campagne particulière (taille des lots définis), l'équipement est décrit par classe (caractéristique telle que la capacité) mais aucun appareil particulier n'est encore désigné. Elle est utilisé au niveau ordonnancement.

Recette de site adaptée à une campagne particulière (taille des lots définis), l'équipement est décrit par classe (caractéristique telle que la capacité) mais aucun appareil particulier n'est encore désigné. Elle est utilisé au niveau ordonnancement.

Recette appliquée à un lot particulier pendant tout son processus de fabrication. Elle correspond à l'exécution temps-réel de la recette principale avec une affectation effective des équipements.

Recette appliquée à un lot particulier pendant tout son processus de fabrication. Elle correspond à l'exécution temps-réel de la recette principale avec une affectation effective des équipements.

Recette de site adaptée à une campagne particulière (taille des lots définis), l'équipement est décrit par classe (caractéristique telle que la capacité) mais aucun appareil particulier n'est encore désigné. Elle est utilisé au niveau ordonnancement.

Recette de site adaptée à une campagne particulière (taille des lots définis), l'équipement est décrit par classe (caractéristique telle que la capacité) mais aucun appareil particulier n'est encore désigné. Elle est utilisé au niveau ordonnancement.

Recette appliquée à un lot particulier pendant tout son processus de fabrication. Elle correspond à l'exécution temps-réel de la recette principale avec une affectation effective des équipements.

Recette appliquée à un lot particulier pendant tout son processus de fabrication. Elle correspond à l'exécution temps-réel de la recette principale avec une affectation effective des équipements.

GRAFCET Niveaux de procédure de la recette Niveaux de procédure de la recette

  • une phase est composée d'une succession de pas permettant de réaliser une fonction élémentaire. Elle est généralement délimitée par des frontières définies par les instants où le déroulement des pas peut être interrompu en toute sécurité pour l'installation et sans dommage pour le produit.• un pas désigne une séquence temporelle ou événementielle d'instructions à exécuter dans un module d'équipement en vue d'accomplir une tâche spécifique. • enfin, une instruction définit une commande élémentaire orientée équipement adressée à un élément ou à un composant afin qu'il effectue un traitement donné.

	est constituée de est constituée de
	3 opérations 3 opérations
	est constituée de 4 phases est constituée de 5 pas est constituée de 4 instructions correspond est constituée de 4 phases est constituée de 5 pas est constituée de 4 instructions correspond	Procédure : Fabrication du PVC Polymérisation du MVC Récupération MVC résiduel Séchage du PVC Procédure : Fabrication du PVC Polymérisation du MVC Récupération MVC résiduel Séchage du PVC Opération : Polymérisation du MVC Préparer le réacteur, Charger le réacteur, Provoquer la réaction, Stopper la réaction. Opération : Polymérisation du MVC Préparer le réacteur, Charger le réacteur, Provoquer la réaction, Stopper la réaction. Phase : charger le réacteur Transférer 50 m 3 d'eau, Ajouter 0,05 m 3 de colloïdes Agiter Ajouter 50 m 3 de MVC Ajouter l'initiateur Phase : charger le réacteur Transférer 50 m 3 d'eau, Ajouter 0,05 m 3 de colloïdes Agiter Ajouter 50 m 3 de MVC Ajouter l'initiateur Polymérisation du MVC Récupération MVC résiduel Séchage du PVC Procédure : Fabrication du PVC Polymérisation du MVC Récupération MVC résiduel Séchage du PVC Opération : Polymérisation du MVC Préparer le réacteur, Charger le réacteur, Provoquer la réaction, Stopper la réaction. Opération : Polymérisation du MVC Préparer le réacteur, Charger le réacteur, Provoquer la réaction, Stopper la réaction. Phase : charger le réacteur Transférer 50 m 3 d'eau, Ajouter 0,05 m 3 de colloïdes Agiter Ajouter 50 m 3 de MVC Ajouter l'initiateur Phase : charger le réacteur Transférer 50 m 3 d'eau, Ajouter 0,05 m 3 de colloïdes Agiter Ajouter 50 m 3 de MVC Ajouter l'initiateur GRAFCET Procédure : Fabrication du PVC GRAFCET GRAFCET
	au signal au signal
		Instruction : Ouvrir V eau signal OUV_Veau Instruction : Ouvrir V eau signal OUV_Veau Instruction : Ouvrir V eau Instruction : Ouvrir V eau signal OUV_Veau signal OUV_Veau

Niveaux de procédure de la recette Procédure : Fabrication du PVC Pas : Transférer 50 m 3 d'eau Ouvrir V eau , Démarrer pompe P 1 Arrêter P 1 si niveau atteint Fermer V eau , Pas : Transférer 50 m 3 d'eau Ouvrir V eau , Démarrer pompe P 1 Arrêter P 1 si niveau atteint Fermer V eau , Procédure : Fabrication du PVC Pas : Transférer 50 m 3 d'eau Ouvrir V eau , Démarrer pompe P 1 Arrêter P 1 si niveau atteint Fermer V eau , Pas : Transférer 50 m 3 d'eau Ouvrir V eau , Démarrer pompe P 1 Arrêter P 1 si niveau atteint Fermer V eau ,

1. Fabrication du produit P selon la réaction (R) dans le réacteur BR1 ou BR2,

	R2 R2
	R1 R1	SYNTHESE 1 SYNTHESE 1	6,5% R1 6,5% R1 6,5% R2 6,5% R2 83% P 83% P	98 % P 98 % P 1% R1 1% R1 1% R2 1% R2	Bidons Bidons de de produit produit
			SEPARATION SEPARATION	CONDITIONNEMENT CONDITIONNEMENT
		SYNTHESE 2 SYNTHESE 2	40% R1 40% R1 40% R2 40% R2
			20% P 20% P
	Cette opération comporte les phases suivantes :
	o charger le réacteur avec n/2 moles de produit R1,
	o préchauffer à 383 K,
	o charger le réacteur de n/2 moles de produit R2 tout en
	régulant la température à 383 K.
	o laisser réagir jusqu'à ce que la composition en produit P atteigne la valeur de 0,8.	Le produit P, stocké dans le réservoir ST4 est ensuite conditionné dans des bidons de 10 litres
	o vidanger tout ou partie du mélange réactionnel.

ρs,i : proportion de s consommée (ρs,i <0) ou produite (ρs,i >0) par i αi,j: partie fixe de la durée d'une tâche i sur une unité j βi,j : partie variable de la durée d'une tâche i sur une unité j H : horizon de temps Demands : demande en s qui doit être satisfaite en fin de campagne

Variables : w(i,j,n) : variable binaire qui associe le début de i sur j à l'évènement n ts(i,j,n) : date début de i sur j à l'évènement n b(i,j,n) : taille du lot de i sur j à l'évènement n ST(s,n) : quantité de s à l'évènement n MS : makespan

.10)

Vers la supervision des procédés

Pour des raisons de sécurité et de performance, les fonctions de surveillance et de supervision ont un rôle majeur en conduite de procédé. Cependant, certains systèmes industriels rendent complexes ces fonctions en raison du nombre croissant de variables process mises en jeu et de la complexité de leurs relations. Le comportement du procédé est alors difficilement interprétable. Pour cette raison, les travaux de recherche en cours visent à concevoir un système de supervision couplé au simulateur dynamique hybride PrODHyS. L'architecture de ce système est montrée sur la figure 22.