David Cantor
email: david.cantor@umontpellier.fr

Émilien Azéma

Philippe Sornay

Farhang Radjai

E Azéma

• F Radjai

Three-dimensional bonded-cell

Keywords: Bonded-cell model, Fragmentation, Discrete element method, Contact dynamics method, Voronoi cell, Weibull statistics

published or not. The documents may come

Introduction

The fragmentation of particles is an important factor in everyday use of granular materials and in processing industries dealing with powders and grains. It modifies the grain size distribution [START_REF] Bandini | The influence of particle breakage on the location of the critical state line of sands[END_REF][START_REF] Fukumoto | Particle breakage characteristics of granular soils[END_REF][START_REF] Mcdowell | The fractal crushing of granular materials[END_REF], packing fraction [15, [START_REF] Hégron | Compaction of a bed of fragmentable UO2 particles and associated acoustic emission[END_REF][START_REF] Lade | Significance of particle crushing in granular materials[END_REF], shear strength [START_REF] Lobo-Guerrero | Discrete element method evaluation of granular crushing under direct shear test conditions[END_REF], yielding surface [START_REF] Miura | Stress-strain characteristics of sand in a particle-crushing region[END_REF][START_REF] Russell | Crushing of particles in idealised granular assemblies[END_REF], and microstructure of granular materials [START_REF] Ma | Modeling the particle breakage of rockfill materials with the cohesive crack model[END_REF]. However, there is presently no clear quantitative understanding of such effects in the rheology and processing of granular materials and, conversely, the properties of the fragments such as their size, shape, and strengths cannot be predicted from the process [START_REF] Guimaraes | Aggregate production: Fines generation during rock crushing[END_REF].

The effects of particle crushing can be investigated experimentally but its origins are much more difficult to trace back to the particle scale. In a similar vein, some continuum mechanics models are able to reproduce phenomenologically the effects of particle fragmentation or degradation [START_REF] Kh | Numerical simulation of particle breakage of angular particles using combined dem and fem[END_REF][START_REF] Cecconi | A constitutive model for granular materials with grain crushing and its application to a pyroclastic soil[END_REF]. Such models are, however, mostly based on idealized description of grain failure. For these reasons, numerical simulations based on the discrete element method (DEM) have been extensively used for a better understanding of the fragmentation process and its impact on the macroscopic behavior. The first DEM simulations of crushable particles were proposed in the late 80s [START_REF] Zubelewicz | Interface element modeling of fracture in aggregate composites[END_REF] and 90s [START_REF] Åström | Fragmentation of grains in a twodimensional packing[END_REF][START_REF] Kun | A study of fragmentation processes using a discrete element method[END_REF][START_REF] Tsoungui | Numerical model of crushing of grains inside two-dimensional granular materials[END_REF]. Initially, these models were two-dimensional and involved simplified (circular) shapes of fragments and their interactions (frictional con-

BCM with contact dynamics 2.1 Voronoi meshing

We perform a three-dimensional meshing of polyhedral particles with the help of the package NEPER [43] and using a centroidal Voronoi tessellation [START_REF] Du | Centroidal Voronoi tessellations: applications and algorithms[END_REF]. In general, a Voronoi tessellation is a division of a space X in R N via a set of points or generators {P j }. The main property of this division is that any point x belongs to a Voronoi cell V i if a distance function d(x, P i ) is minimal between all the generators [START_REF] Okabe | Spatial tessellations: concepts and applications of voronoi diagrams[END_REF]. A compact definition of a Voronoi cell may be written as follows:

V i = {x ∈ X | d(x, P i ) < d(x, P j ) ∀ j = i}.
(

) 1
In three dimensions, the Voronoi tessellation produces always convex polyhedra sharing common faces if the region X is also convex. Furthermore, a centroidal Voronoi tessellation has the property that the generators are also the centroids of the cells. The points {P j } are usually located randomly and, as a consequence, the locations of those points differ from those of the corresponding cell centroids {c j }. One possible technique to reduce the distances d(P j , c j ) is known as the Lloyd's method. This iterative procedure may be described as follows:

1. Define the space X, 2. select points {P j }, 3. define the corresponding cells {V j }, 4. compute the centroid of each cell {c j }, 5. replace points {P j } with points {c j } and recompute {V j }, 6. repeat from step 4 until a convergence criterion is satisfied.

The convergence criterion that is used in this work may be written as follows:

n cl j=1 d(P j , c j ) n cl j=1 r j ≤ λ, ( 2
)
where n cl is the number of cells, and r j is the average radius of cell j. This ratio should be below or equal the tolerance λ at which we stop the iterations. The quantity λ is actually a quite interesting parameter because, as it is shown below, it may be interpreted as a measure of disorder of the cell shapes and their distributions. In general, values of λ 0 produce less elongated shapes and tacts). 2D simulations are now capable of reproducing the evolution of cracks and complex geometries of fragments [START_REF] Nguyen | Evolution of particle size distributions in crushable granular materials[END_REF][START_REF] Nguyen | Bonded-cell model for particle fracture[END_REF].

With increasing computational power, three-dimensional models have also been introduced recently. These models can be classified at least in two main groups. The first group, which we call Replacing Method, consists in substituting each broken particle when typically a load threshold is reached by a set of smaller particles [START_REF] Ciantia | An approach to enhance efficiency of dem modelling of soils with crushable grains[END_REF][START_REF] Zhou | Macromicro responses of crushable granular materials in simulated true triaxial tests[END_REF]. However, a nonconservation of mass and volume may occur as the replacements are made and artificial overlapping between particles is sometimes allowed, producing local elastic overstresses.

The second group, which was coined Bonded-Cell Method (BCM) in Ref. [START_REF] Nguyen | Bonded-cell model for particle fracture[END_REF], is based on the construction of particles by agglomerating a set of potential fragments. These small fragments are held together by diverse mechanisms that assign a bonding strength at each contact. The bonding of spheres is the most wide-spread model to simulate crushing in 3D. However, the spheres are obviously unable to reproduce the broad shape variability of fragments. Moreover, the apparent volume of a sphere agglomerate is considerably larger than those of its fragments. In order to circumvent these problems, some authors have used polygonal or polyhedral subparticles or cells generated by Voronoi tessellation [START_REF] Galindo-Torres | Breaking processes in three-dimensional bonded granular materials with general shapes[END_REF]. These cells pave the whole volume of the particle so that the volume is conserved during particle fracture and fragmentation. Nevertheless, intercell contacts are modeled by a linear force law as that between spherical subparticles. This is clearly an unphysical approximation since the contacts extend along a surface between cells and a realistic model must account for the fact that the treatment of such contacts need at least three displacement variables at three independent contact points.

In this paper, we introduce a three-dimensional BCM approach in the framework of the contact dynamics (CD) method. In 3D, each particle is composed of cells of irregular polyhedral shape. The framework of the CD method has the advantage of accounting for the correct kinematics of faceface contacts between fragments [START_REF] Jean | The non-smooth contact dynamics method[END_REF][START_REF] Radjai | Contact dynamics as a nonsmooth discrete element method[END_REF]. Furthermore, the treatment of cells as rigid particles leads to numerically efficient simulations in 3D. We also introduce a general cohesive law governing cell interactions before fracture. Using this approach, we investigate the fracture properties of a single particle subjected to diametral compression.

This paper is organized as follows: In Sect. 2.3,we present the 3D BCM as well as the treatment of intercell contacts within the CD method. In Sect. 3, numerical simulations of diametral compaction are introduced. In Sect. 3.2, we investigate the effects of different geometrical and mechanical parameters on the particle strength. Finally, in Sect. 4, a brief summary and perspectives of this work are presented. By construction, the sizes of the initial aggregates and cells determine the upper and lower bounds of the fragment size distributions during crushing simulations. In this sense, for DEM simulations it is important to optimize the number of cells per aggregate in order to ensure the representativity of the fragments sizes for a convenient computational time.

Cell interactions

The contact between two polyhedral shapes can have three different geometries: contact point, contact line, and contact surface, as shown schematically in Fig. 3. The determination of the contacts between two polyhedral cells is performed in several steps. First, a "bounding box" method is used to compute a list of neighboring pairs. Then, for each pair, the overlaps are calculated through the so-called "common plane" method introduced by Cundall [17]. This is an iterative method based on the perturbation of the orientation of the normal vector. The process is initialized by a first guess based on the vector joining the centers of the two polyhedra. This detection procedure is fairly rapid when correctly optimized [37,[START_REF] Nezami | Shortest link method for contact detection in discrete element method[END_REF]. The interactions may include The contact line, such as that produced by a edge-face interaction, can be represented by two distinct points, whereas a surface (face-face) contact can be replaced by three points since they involve an equivalent number of geometrical unilateral constraints between two faces. Note that vertex-vertex or vertex-edge contacts are rare but, when they occur, the common plane method is able to provide a normal direction based on geometrical criteria.

In the Bonded-Cell Method, we assume that rigid cells interact only via cohesive forces through face-face contacts. The normal adhesion threshold f c n depends on the area S of the contact. Since the contact plane is represented by three points, the tensile threshold for each contact point is given by f c n = C n S/3, where C n is the internal cohesion of the material. In the same way, the shear threshold for each contact point in the contact plane is given by f c t = C t S/3, where C t is the contact shear strength. C n and C t are assumed to be independent strength parameters of the material. Once either of the critical stresses is reached, the intercell bond breaks irreversibly and the interaction turns into a frictional contact with coefficient of friction μ cl . Note that, since the cells are rigid, a cohesive face-face contact can not kinematically break until the three points representing the interface do so, i.e., when either of the two thresholds is reached for each of the three points.

The normal force f n and normal velocity u n at each contact point satisfy the following relations [START_REF] Radjai | Contact dynamics as a nonsmooth discrete element method[END_REF]

u n > 0 ⇒ f n = 0 u n = 0 ⇒ f n + C n S/3 > 0 (3)
which simply expresses the unilateral nature of the contact between the two cells at their contact points in the presence of a cohesive force. In the same way, along the tangential direction, we have the following relation between the tangential force f t and the tangential velocity u t :

⎧ ⎨ ⎩ u t > 0 ⇒ f t = -C t S/3 u t = 0 ⇒ -C t S/3 ≤ f t ≤ C t S/3 u t < 0 ⇒ f t = C t S/3 . ( 4
)
These "complementarity relations" are represented by two graphs shown in Fig. 4a,b.

The above relations should be taken into account together with the equations of motion in order to determine both the velocities and forces at all contacts in the system. We use the contact dynamics (CD) method, which, as briefly discussed below, allows for exact prescription of these relations based on an implicit time-stepping integration scheme. To use molecular dynamics (MD) or common DEM, it is necessary to regularize the above relations so that the classical explicit integration schemes can be applied.

At each time step, the "critical" contact points, where u n > 0 (separation) or |u t | > 0 (sliding), break irreversibly and the contact point turns into a cohesionless frictional contact governed by the following relations:

(5)

⎧ ⎨ ⎩ u t > 0 ⇒ f t = -μ cl f n u t = 0 ⇒ -μ cl f n ≤ f t ≤ μ cl f n u t < 0 ⇒ f t = μ cl f n . ( 6
)
The graphs of these relations are shown in Fig. 4c,d.

Contact dynamics

The above complementarity relations were prescribed within the CD method [START_REF] Jean | The non-smooth contact dynamics method[END_REF][START_REF] Moreau | Some numerical methods in multibody dynamics: application to granular[END_REF][START_REF] Radjai | Contact dynamics as a nonsmooth discrete element method[END_REF][START_REF] Moreau | Numerical investigation of shear zones in granular materials[END_REF]. In this method, the rigid-body equations of motion are integrated for all particles, which are treated as undeformable bodies, by taking into account the kinematic constraints resulting from the complementarity relations. An iterative algorithm similar to Gauss-Seidel iterations is used to determine the contact forces and particle velocities simultaneously at all potential contacts. Since this algorithm is implicit, the contact velocities u n and u t entering the complementarity relations can represent the velocities u + n and u + t at the end of a time step. This is equivalent to the assumption that the coefficients e n and e t of normal and tangential restitution are both zero. However, it is also possible to use the complementarity relations with the following weighted means:

u n = e n u - n + u + n 1 + e n ( 7
)
u t = e t u - t + u + t 1 + e t , ( 8
)
where u - n and u - t are the normal and tangential contact velocities at the beginning of a time step. With complementary relations involving these weighted velocities, it can be shown that in a binary collision between two particles the postcollisional velocities are indeed given by u + n = -e n u - n and u + t = -e t u - t . In this way the contact laws involve five independent parameters: coefficient of friction μ cl , coefficient of normal restitution e n , coefficient of tangential restitution e t , normal cohesive strength C n , and tangential cohesive strength C t .

An implicit time-stepping scheme makes the CD method unconditionally stable. In contrast to the molecular dynamics method, in the CD method the tiny numerical overlaps between particles are used for contact detection but they do not represent an elastic deflection. For this reason, the time step can be larger than that in the MD method or DEM. The CD method has been extensively employed for the simulation of granular materials in 2D and 3D [START_REF] Azéma | Stress-strain behavior and geometrical properties of packings of elongated particles[END_REF][START_REF] Azéma | Force chains and contact network topology in sheared packings of elongated particles[END_REF][START_REF] Azéma | Force transmission in a packing of pentagonal particles[END_REF][START_REF] Azéma | Rheology of 3D packings of aggregates: microstructure and effects of nonconvexity[END_REF][START_REF] Azéma | Quasistatic rheology, force transmission and fabric properties of a packing of irregular polyhedral particles[END_REF][START_REF] Bratberg | Dynamic rearrangements and packing regimes in randomly deposited two-dimensional granular beds[END_REF][START_REF] Estrada | Identification of rolling resistance as a shape parameter in sheared granular media[END_REF][START_REF] Nezamabadi | Implicit frictional-contact model for soft particle systems[END_REF][START_REF] Nouguier | Force fluctuations in a pushed granular material[END_REF][START_REF] Quezada | Stability, deformation, and variability of granular fills composed of polyhedral particles[END_REF][START_REF] Radjai | Force Distributions in Dense Two-Dimensional Granular Systems[END_REF]48,[START_REF] Ries | Shear zones in granular media: Three-dimensional contact dynamics simulation[END_REF][START_REF] Saussine | Modelling ballast behaviour under dynamic loading. Part 1: A 2D polygonal discrete element method approach[END_REF][START_REF] Saint-Cyr | Rheology of granular materials composed of nonconvex particles[END_REF][START_REF] Staron | Multi-scale analysis of the stress state in a granular slope in transition to failure[END_REF][START_REF] Staron | Preavalanche instabilities in a granular pile[END_REF][START_REF] Taboada | Rheology, force transmission, and shear instabilities in frictional granular media from biaxial numerical tests using the contact dynamics method[END_REF][START_REF] Topin | Collapse dynamics and runout of dense granular materials in a fluid[END_REF]. We implemented the Bonded-Cell method in the software LMGC90 [START_REF] Dubois | LMGC90 wiki page[END_REF]. Video samples of diametral compression tests can be found at www.cgp-gateway.org/ref035.

u n > 0 ⇒ f n = 0 u n = 0 ⇒ f n > 0 and

Diametral compression of an icosahedron

We consider icosahedral particles of mean radius R compressed between two platens, as illustrated in Fig. 5. A gradually increasing vertical force f a is applied on the top platen. The platens directly touch two opposite hexagonal faces of the icosahedron. Hence, the force f a is always exerted on this supporting hexagonal surface belonging to the periphery of the particle. This surface may be shared between several cells depending on the number n cl of cells inside the icosahedron. As f a is increased from zero, the contact forces inside the particle increase proportionally to f a but, since the particle is rigid, no fracture occurs until a threshold f c is reached. At this point, an unstable failure of the particle occurs and the particle breaks into several fragments. Several snapshots of a particle at incipient crushing are shown in Fig. 6. The subsequent fragmentation is a dynamic process and can be investigated under controlled displacement of the top platen. Here we are interested in the initiation of failure and the breaking force f c , which represents the particle strength.

The intercell breaking strengths C n and C t having the dimensions of a stress, the particle strength should be expressed as a stress threshold σ c . However, the stress field inside a particle subjected to diametral compression is nonuniform [START_REF] Chau | Spherically isotropic elastic spheres subject to diametral point load strength test[END_REF]60]. The average vertical stress in the largest section of the particle π R 2 is f a /π R 2 . The highest tensile stress inside the particle can be much larger than this stress. In our meshed particles, even with a large number of cells, the disordered configuration of the cells results in a strongly inhomogeneous distribution of intercell forces. Figure 7 shows an example of intercell forces in a particle composed of 150 cells. We observe both compressive forces along the vertical direction and tensile forces along the horizontal direction. We can compute the average vertical stress from the forces by applying the following well-known micromechanical expression:

σ a = 1 V α∈V f α z α z , ( 9
)
where V is the volume of the particle, α denotes the contact labels between the cells, f z is the vertical component of the intercell force, and z is the vertical component of the branch vector joining the centers of cells. This average stress is proportional to f a /π R 2 by a factor of 1.2 and it reaches its maximum value for f c . Note that, to avoid inhomogeneities induced by gravity, the latter is set to zero although its value is by orders of magnitude below the breaking threshold f c .

Effect of intercell strength parameters

Since the stresses and forces inside the particle scale with the applied force, it is generally expected that the particle strength is a linear function of the tensile strength C n of intercell contacts. However, we expect the shear strength C t to play a crucial role since an intercell contact fails at this threshold (when sliding occurs) and both C n and C t vanish upon failure. Generally, the shear strength C t is a function of C n depending on the material [START_REF] Barton | The shear strength of rock and rock joints[END_REF][START_REF] Barton | Shear strength criteria for rock, rock joints, rockfill and rock masses: Problems and some solutions[END_REF][START_REF] Stoller | On the relationship between uniaxial yield strength and resolved shear stress in polycrystalline materials[END_REF]. However, in view of exploring the respective roles of C t and C n , we varied these parameters independently in the range between 0.25 and 8 MPa.

Figure 10 shows the evolution of the particle strength σ c as a function of C n for different values of C t . Each data point is an average over 10 independent simulations. We observe

Effects of meshing on the particle strength

We performed a series of simulations with a particle of constant size but varying number n cl of cells and degree of disorder λ. In these simulations, the intercell strengths are kept constant and equal C n = C t . The intercell coefficient of friction, which becomes effective only when an intercell contact fails, is set to μ cl = 0.4. Figure 8 Despite the large variability of the data, specially at small numbers of cells, we see that the particle strength increases on the average with n cl and it becomes independent of it beyond 75 cells. We also observe that this steady value of the particle strength is nearly 1.25 times the intercell strength C n .The high variability of the particle strength with low number of cells reflects the very high inhomogeneity of forces for such a low number of intercell contacts. The particle may break just as a result of the failure of a single intercell contact. For this reason, the average values in this range for 10 independent simulations do not represent statistically well-defined averages.

In order to evaluate the effect of the degree of disorder in cell configurations, we used particles composed of 150 cells and an increasing degree of disorder λ. Figure 9 displays the normalized strength as a function of λ. We see that, as λ increases from zero, the particle strength rapidly declines. The particle strength is nearly constant beyond λ = 0.2. The transition between these two limits seems to be nonlinear. This variation can be explained by remarking that ordered cell structures are formed for λ close to zero where the cells The solid lines are power-law fits to the data a nonlinear dependence of the particle strength on both cell strength parameters. At low values of C n , σ c is rather low but it varies by a factor of 20 as C t increases to its largest value. At largest values of C n , on the other hand, σ c can take values as much as two times larger but varies by only a factor of 4 as C t increases in the same range. This shows that the effect of C t directly depends on the value of C n . In the same figure, we also see that the particle strength for C n = C t increases nearly linearly with the latter with a proportionality coefficient of 1.25.

The nonlinear trends observed in Fig. 10 are clearly a consequence of the variation of failure mode in response to the relative values of C t and C n . In particular, we expect that at low values of C t with respect to C n , particle breakage by intercell slip requires less energy dissipation and is thus more likely to occur, whereas at high values of C t with regard to C n , failure by tensile rupture should prevail. This picture is consistent with the evolution of the normalized particle strength σ c /C n as a function of the ratio C t /C n plotted in Fig. 11. First, we see that all data points from simulations nicely collapse on a single curve. Secondly, the data are well fitted by two straight lines of slopes α and β representing power laws with exponents α and β for C t < C n and C t > C n , respectively, over nearly three decades. Hence, to a good approximation, we have

σ c C n = 1.25 × ⎧ ⎨ ⎩ C t C n α C t ≤ C n C t C n β C t > C n , ( 10
)
where α 0.5 and β 0.8. Note how these two relations are reduced to a linear function of C n for C t /C n = 1 with the same prefactor 1.25 as in Fig. 10.

Equation [START_REF] Barton | The shear strength of rock and rock joints[END_REF] with the values of the exponents α and β can be written as σ c = 1.25 √ C n C t in the range C t ≤ C n . This form is a generalization of the observed linear relation in the case C n = C t . In the range C t > C n , we have σ c

1.25C 0.2 n C 0.8 t . This relation implies that the particle strength increases faster here with C t than with C n . The exponent β is so close to 1 that, in practice, σ c varies almost linearly in response to C t . This seems to contrast the expected picture that in this regime, where C t is above C n , tensile debonding is more likely to occur, implying thus that the particle strength should be mainly controlled by C n . But the larger values of the intercell shear strength compared to tensile strength is a special condition in which the number of possible configurations accommodating cell motions, and hence failure, declines as C t increases. In other words, tensile forces are activated to a lower extent since force balance over cells is increasingly ensured by shear forces. Most of time, a mechanically compatible displacement field of the cells at incipient failure involves both tensile debonding and slip at a subset of intercell contacts. This is very similar to the effect of friction in granular materials. Extensive simulations indicate that in packings of polyhedral particles, the friction force mobilization and its anisotropy play a much more important role for the shear strength than in a packing of spherical particles [START_REF] Azéma | Quasistatic rheology, force transmission and fabric properties of a packing of irregular polyhedral particles[END_REF]. This effect is reminiscent of "frictional interlocking" in analogy to "geometrical interlocking" of strongly nonconvex particles [START_REF] Jaeger | Celebrating soft matter's 10th anniversary: Toward jamming by design[END_REF]. The values of the exponents in the above two regimes reflect thus the discrete texture of the particle, steric exclusions between rigid cells and frictional interlocking.

Strength variability

The foregoing analysis of particle strength as a function of intercell strength refers to the measured average values of the particle strength for several independent cell configurations. However, particle-to-particle fluctuations occur to a large extent, and the particle can survive under compressive stresses above its mean strength. As in rocks and agglomerates of fine particles, such fluctuations reflect the presence and distribution of defects and stress concentration factors, which vary with particle size and/or the size of defects compared to particle size [START_REF] Mcdowell | On the micromechanics of crushable aggregates[END_REF][START_REF] Nguyen | Bonded-cell model for particle fracture[END_REF]. In our simulations, the particle size is kept constant and the number n cl of cells is varied in a limited range such that our largest spheres contain only 250 cells. This is too small to allow for a systematic analysis of size effect as the particle size increases (except for particles containing below 75 cells and showing slightly lower strength as discussed in Sect. 3.2). For this reason, we focus here on strength variability only for particles composed of 150 cells.

Figure 12 shows the survival probability P s (the probability that the particle does not break under a load σ a ) of normalized particle stress σ a /C n in the case C n = C t . The log-log plot is consistent with Weibull probability distribution:

P s (σ a /C n ) = e -(σ a /σ 0 ) m , (11)
to obey Tresca's failure criterion. Specifically, the bonding strength is controlled by critical stresses along normal and shear directions between cells. The cells were generated by means of a three-dimensional Voronoi meshing in which the control parameters are the number of cells and their degree of disorder.

A systematic study was undertaken to analyze the effect of various parameters on the particle strength defined as the compressive strength of the particle under diametral loading between two platens. For this purpose, particles of icosahedral shape were loaded between two platens. We showed that particle strength increases with the number of cells and levels off at around 1.25 times the tensile intercell strength. On the other hand, as disorder increases, the particle strength rapidly declines to a well-defined value. The particle strength in those asymptotic states was shown to define two distinct regimes as a function of the ratio of intercell shear strength to tensile strength. In both regimes, the particle strength varies as a power law of both shear and normal intercell strengths. We argued that this nonlinear behavior is a consequence of the evolution of failure modes at the cell scale as the above ratio increases. In particular, the particle strength depends not only on the intercell tensile strength but also crucially on the intercell shear strength, which prevails when the strength ratio is above 1, and amplifies interlocking effect between cells. Finally, we also showed that particle strength fluctuations are consistent with Weibull statistics.

Our bonded-cell model in 3D and its behavior analyzed here under diametral loading provide the first applications of this model in the framework of the contact dynamics method with particles and cells treated as perfectly rigid bodies. This model is simple and capable of reproducing the angular shapes of fragments while conserving the total volume. Further analyses are underway to highlight the initiation of particle failure at the cell scale, the effects of cell disorder on strength variability, and finite size effect as the number of cells is increased. This method can be applied to simulate large samples of crushable particles in 3D.

where σ 0 = kC n is a scaling stress parameter and m is the Weibull modulus. Here, we have m 2, which corresponds to a Gaussian distribution around the mean strength σ c , and k 1.4. This is below experimental values, which typically vary in the range from 3 to 10 [START_REF] Hardin | Crushing of soil particles[END_REF][START_REF] Lade | Significance of particle crushing in granular materials[END_REF][START_REF] Mcdowell | On the micromechanics of crushable aggregates[END_REF]. But m is expected to increase with disorder parameter λ and particle size [START_REF] Nguyen | Bonded-cell model for particle fracture[END_REF]. This value of m is an interesting observation since it shows that normal fluctuations should be expected in this range of the rather small numbers of cells.

Summary and perspectives

In this paper, a three-dimensional model for grain crushing was proposed in the framework of the Contact Dynamics Method. This model belongs to the class of bonding methods as each particle is discretized into a collection of polyhedral cells. The bonding strength between fragments was assumed

Figures 1 and 2

 2 display several examples of meshing of a truncated icosahedron, a semiregular polyhedron with 60 vertices and 21 faces, for λ = 0.25 and two different numbers of cells, and for n cl = 250 and two values of λ, respectively.

Fig. 1

 1 Fig. 1 Examples of icosahedral particles for λ = 0.25 and meshed by a 5 and b 250 cells

Fig. 2

 2 Fig. 2 Examples of icosahedral particles composed of 250 Voronoi cells and disorder parameters λ = 0.005 (a) and λ = 0.5 (b)

Fig. 3

 3 Fig. 3 Examples of contact types that may occur between polyhedral cells: a Point, b line, and c surface

Fig. 4

 4 Fig. 4 Complementarity relations along the normal and tangential directions of the contact frame for cohesive contacts (a, b) and broken contacts (c, d) respectively

Fig. 5

 5 Fig. 5 Schematic representation of diametral compression test on a particle of average radius R

Fig. 6

 6 Fig. 6 Snapshots of a particle composed of 75 cells undergoing incipient breakage under the action of a vertical force

Fig. 8 9

 8 Fig. 8 Evolution of particle strength with the number of cells. Each point represents an average value over 10 independent cell configurations. The error bars represent minimum and maximum values

Fig. 7

 7 Fig. 7 Intercell forces in a particle composed of 150 cells subjected to diametral compression. Line thickness is proportional to the force. The compressive and tensile forces are in red and blue, respectively

  shows the vertical strength σ c , normalized by normal strength C n , as a function of n cl . Each point represents the mean value for 10 simulations performed with particles composed of cells with the same degree of disorder λ = 0.25 and the same strength, but independent Voronoi tessellations for the generation of the cells. The maximum and minimum values are shown by error bars.

Fig. 10 Fig. 11

 1011 Fig. 10 Evolution of particle strength as a function of intercell tensile strength C n for different values of intercell shear strength C t . All stresses are in MPa. The circles indicate the data points for C n = C t

Fig. 12

 12 Fig.[START_REF] Bratberg | Dynamic rearrangements and packing regimes in randomly deposited two-dimensional granular beds[END_REF] Survival probability of particles as function of the mean applied stress for particles composed of 150 cells

Acknowledgments

This work was financially supported by a research grant awarded by the French Alternative Energies and Atomic Energy Commission (CEA). Farhang Radjai would also like to acknowledge the support of the ICoME2 Labex (ANR-11-LABX-0053) and the A*MIDEX projects (ANR-11-IDEX-0001-02) cofunded by the French program Investissements d'Avenir, managed by the French National Research Agency (ANR).


We present detailed investigations on single crystals of quaternary EuRhAl 4 Si 2 and EuIrAl 4 Si 2 . The two compounds order antiferromagnetically at T N1 = 11.7 and 14.7 K, respectively, each undergoing two magnetic transitions. The magnetic properties in the ordered state present a large anisotropy despite Eu 2+ being an S-state ion for which the single-ion anisotropy is expected to be weak. Two features in the magnetization measured along the c-axis are prominent. At 1.8 K, a ferromagnetic-like jump occurs at very low field to a value one third of the saturation magnetization (1/3 M 0 ) followed by a wide plateau up to 2 T for Rh and 4 T for Ir-compound. At this field value, a sharp hysteretic spinflop transition occurs to a fully saturated state (M 0 ). Surprisingly, the magnetization does not return to origin when the field is reduced to zero in the return cycle, as expected in an antiferromagnet. Instead, a remnant magnetization 1/3 M 0 is observed and the magnetic loop around the origin shows hysteresis. This suggests that the zero field magnetic structure has a ferromagnetic component, and we present a model with up to third neighbor exchange and dipolar interaction which reproduces the magnetization curves and hints to an "up-up-down" magnetic structure in zero field.

The isothermal magnetization in an antiferromagnet at temperature T < T N , where T N is the Néel temperature, increases with field from zero to full saturation at the spin-flip transition, when all the moments are aligned parallel to applied field. In the particular case where the field is parallel to the easy-axis along which the spins are pinned, the magnetization in an ideal collinear bipartite antiferromagnet is zero (T < < T N ) until the field induces a spin orientation, or spin-flop/metamagetic transition, leading to a sharp increase/saturation in magnetization. Often the spin-flop/metamagnetic transition is first order in nature as revealed by magnetic hysteresis when the polarity of the field is reversed. In all cases that we are aware of, the magnetization of an antiferromagnet decreases to zero when the field is reduced to zero. In contrast, the magnetization of a ferromagnet shows a hysteretic behaviour around the origin, characterized by a finite remanence and a coercive field. Among other factors, the single-ion magnetocrystalline anisotropy plays an important role in the coercive field of a ferromagnet. Here, we report, to our knowledge the very first case, the observation of hysteresis around the origin in two antiferromagnets, which are thought to present a peculiar magnetic structure. The unusual nature of our observation is further heightened by the fact that the two antiferromagnets EuTAl 4 Si 2 (T = Rh and Ir) derive their magnetism from Eu 2+ which is an S-state ion (S = 7/2; L = 0). For S-state rare earth ions like Eu 2+ and Gd 3+ , the single-ion magnetocrystalline anisotropy is usually quite small.

The synthesis of single crystals of new quaternary compounds EuTAl 4 Si 2 (T = Rh and Ir), using the Al-Si binary eutectic as flux, has recently been reported [START_REF] Maurya | Crystal and Electronic Structure of the Quaternary Magnetic EuTAl 4 Si 2 (T = Rh and Ir) Compounds[END_REF] . The two compounds were found to adopt an ordered derivative of the ternary KCu 4 S 3 -type tetragonal, tP8, P4/mmm structure. The Eu, Al, Si and T atoms occupy respectively the K(1a), Cu(4i), S2(2h) and S1(1a) Wyckoff sites, which leads to quaternary and truly stoichiometric 1:1:4:2 compounds. The local symmetry at the Eu site is fourfold axial (4/mmm). Preliminary magnetization data [START_REF] Maurya | Crystal and Electronic Structure of the Quaternary Magnetic EuTAl 4 Si 2 (T = Rh and Ir) Compounds[END_REF] , taken with the applied field parallel to [001] direction, revealed a divalent state for Eu ions in both compounds and a magnetic transition at ≈ 12 and ≈ 15 K in the Rh and Ir analogs, respectively. Electronic structure calculations using the local spin density approximation 1 are in conformity with the divalent state of Eu ions.

In the present work, we probe in detail the magnetic behavior of these two intermetallic compounds along the principal crystallographic directions of the tetragonal cell, using the techniques of magnetization, heat capacity and electrical transport. Additional information about the magnetic transitions is derived from 151 Eu Mössbauer spectroscopy. We find a rather unusual and highly anisotropic magnetic response, which points to a probable ferrimagnetic-like spin arrangement, scarcely observed in intermetallics with divalent Eu. We show that such a ground magnetic structure can be predicted using a mean field model with three exchange integrals. ( )

Results

Magnetic susceptibility.

χ μ θ ( ) = - , ( )
T N k T 3 1 B p eff 2
where N is the Avogadro number, almost in the entire paramagnetic region. The effective moment μ eff and the Curie-Weiss temperature θ p are listed in exchange couplings, with opposite signs. This is probably brought about in these two metallic compounds by the oscillatory RKKY (Ruderman-Kittel-Kasuya-Yosida) exchange interaction. An appreciable anisotropy in the susceptibility, plotted between 1.8 and 20 K in the insets of Fig. The susceptibility χ a (T) of EuRhAl 4 Si 2 along the [100] direction (inset, Fig. 1(c)) shows a peak at 12 K followed by a minor kink at 10.4 K (which correspond closely with the peak in zero field heat capacity at 11.7 K and the small kink at 10.4 K, vide infra) with nearly overlapping ZFC and FC data, characteristic of an AF transition. In EuIrAl 4 Si 2 , (inset, Fig. 1(d)), the corresponding data show a peak at nearly 15 K, followed by a barely discernible anomaly in the 13-14 K range (which also corresponds closely with the sharp peak in zero field heat capacity at 14.7 K, followed by a peak at 13.2 K, vide infra). Data taken with H ‖ [110] (not shown) reflect a nearly isotropic magnetic response in the ab-plane in both compounds.

On the other hand, along the tetragonal c-axis [001], the susceptibility χ a (T) exhibits ferromagnetic-like behavior, with its peak value more than an order of magnitude larger than the peak value of χ a (T) (main panel Figs. 1(c,d)). There is a substantial thermomagnetic irreversibility between ZFC and FC runs, which decreases as the applied field is increased (data at higher fields not shown here for the sake of clarity). Phenomenologically, the ZFC/FC behavior is similar to that seen in ferromagnets with a sizeable coercivity. The upturn in χ c (T) begins at nearly the same temperature as the peak in χ a (T). Data were taken in both field-cooled-cooling (FCC) and field-cooled-heating (FCH) modes. The observation of hysteresis in these two plots is in line with the ferromagnetic character of the transition for H ‖ [001].

The behavior of the in-field magnetization at 1.8 K, measured in the ab-plane along [100] and [110], and along [001] (Figs. 2(a,b)), is in conformity with the χ(T) data depicted in Figs. 1(c,d). For H in the ab-plane, the magnetisation in EuRhAl 4 Si 2 increases almost linearly, reaching a (spin-flip) saturation field near 5.5 T, with a saturated Eu 2+ moment M 0 = 7 μ B , as expected for a g = 2, S = 7/2 ion. This behavior is characteristic for an AF structure where the moments are perpendicular to the field. At higher temperatures, both the spin-flip field value and the saturation magnetization decrease (data not shown). By contrast, for H ‖ [001], there is a jump of the magnetization at very low field to a value 2.3 μ B , which is almost exactly 1/3 of M 0 . The magnetization plateau is broken near 1.7-2 T by a hysteretic spin-flop jump to saturation. This behavior tells that the Eu moment probably lies along the c-axis, and the low field jump hints to the presence of an average ferromagnetic component along c in the zero field structure. The 1/3 M 0 value suggests that the latter could consist in a motif of 3 spins along c, two pointing up and one pointing down. The zero field magnetic structure would then consist in ferromagnetic planes with moments along c in an "up-up-down" sequence. The inset of Fig. 2(a) shows that the width of the hysteresis associated with the spin-flop jump decreases and it occurs at lower fields with the increase of temperature. The Ir sibling (Fig. 2(b)) shows exactly the same behavior, but with a spin-flop field of 3.9-4 T along [001] and a spin-flip field of 9.2 T perpendicular to [001]. Figures 2(c,d) show the hysteretic behaviour of the magnetisation in EuRhAl 4 Si 2 and EuIrAl 4 Si 2 respectively for magnetisation jump at H close to origin. Surprisingly, there is a sizeable remanent magnetization of 2.3 μ B /f.u. at 1.8 K and a large coercive field of nearly 0.18 T and 0.25 T in the Rh and Ir analogue, respectively. The width of the loop depends sensitively on the temperature and it has decreased by more than half at 3 K, compared to its value at 1.8 K and by 6 K the hysteresis has practically vanished. We suggest hysteresis is due to the presence of a ferromagnetic component in the zero field structure of these compounds, the possible domains having their mean magnetization pointing either up or down. The coercive field would then correspond to a complete reversal of the domains magnetization along the direction of the field. The primary sources of coercivity are the magnetocrystalline anisotropy and the dipolar field, which has been shown to play an important role in Eu compounds with low transition temperatures 2 . The rather high coercive field of 0.2 T at 1.8 K points thus to an unusually high crystalline anisotropy for an L = 0 divalent Eu ion, for which it vanishes in principle at first order. This is supported by the modelling of the magnetic behavior presented in the discussion section. Magnetic characteristics of EuRhAl 4 Si 2 and EuIrAl 4 Si 2 are summarized in Table 2.

Heat Capacity. The heat capacity of the two compounds is shown in Fig. 3. We have also plotted the heat capacity of the corresponding La analogs, which serves as a measure of the phonon contribution assuming that the phonon spectra in the La and corresponding Eu compound are identical. In EuIrAl along the major crystallographic directions.

an extremely sharp peak with a peak value of nearly 60 J/mol K and peak position of 14.7 K (Fig. 3(a)), indicative of a first order phase transition from the paramagnetic state to an intermediate magnetically ordered state, is followed by a second transition near 13.2 K marked by a small anomaly at that temperature. The heat capacity of EuRhAl 4 Si 2 (Fig. 3(b)) is apparently marked by only one peak at 11.7 K which is relatively less sharp and has a lower peak height. However, a small kink at 10.4 K is discernible; both are marked by arrows in the figure. The 151 Eu Mössbauer spectra, to be described below, also show the presence of two magnetic transitions in both compounds. The magnetic contribution to the heat capacity, C mag , was obtained by subtracting the normalized heat capacity of the corresponding La compound, following the procedure described in Ref. 3. The entropy S mag calculated by integrating C mag /T against the temperature is also shown in Figs. 3(a,b). Nearly 75% of the full entropy R In 8 is recovered at the transition temperature, and the remaining by 40 K.

Electrical Resistivity and Magnetoresistance. Magnetic Field (T) The fall in the resistivity at T N1 , due to the gradual loss of spin disorder scattering below the magnetic transition is clearly visible, while T N2 is marked by only a barely discernible kink in the ρ(T) plot but clearly observed in the derivative plot as shown for EuIrAl 4 Si 2 . To within the precision of our measurements we did not detect any hysteretic behavior of the resistivity when the temperature was cycled up and down across the magnetic transition. The variation of resistivity with temperature in the magnetically ordered state in applied magnetic fields, and the magnetoresistance MR derived from the data can be qualitatively understood on the basis of the magnetization data discussed above. Figures 4(c,d) show the resistivity below 25 K in EuIrAl 4 Si 2 with J ‖ [100] and H ‖ [010] and [001] at selected values of field between 0 and 8 T. In Fig. 4(c) the main anomaly in the resistivity occurring at T N1 shifts towards lower temperatures with increasing field as expected for an antiferromagnet. The anomaly almost disappears at 8 T as the spin-flip field value is approached.

The Rh analog also shows similar behavior except that the anomaly at T N1 vanishes at a lower field consistent with a lower spin-flip field in this compound when the field is applied in the ab-plane. In Fig. 4(d), which shows the data for H ‖ [001], the main anomaly at T N1 again shifts to lower temperature with field. Additionally, the resistivity at 2.75, 3, 3.25, 3.5 and 3.75 T shows an upturn at a temperature that is dependent on the applied field. An upturn, prima-facie suggests a gap-opening; neutron scattering is required to probe if there is a change in the magnetic symmetry with field. The magnetoresistance defined as MR = Δ R/R(0), where Δ R = R(H) -R(0), is shown in Figs. 4(e,f) for EuIrAl 4 Si 2 . The MR of the Rh-analogue is qualitatively similar (data shown in Supplementary information). The MR of an antiferromagnetic material below T N is expected to be positive as the applied field tends to break the antiferromagnetic couplings thereby increasing the spin disorder scattering. Indeed, the MR is positive below T N in both the compounds and increases with increasing field. Near the field-induced spin-flip transition, the MR tends to be field independent (Fig. 4(e)). However, at higher fields the positive MR keeps on increasing though at a somewhat slower rate. In this context it may be mentioned that the MR of non-magnetic reference compound LaIrAl 4 Si 2 is positive and nearly 20% at 14 T; indicating a positive MR due to background conduction electrons. With the increase of temperature the spin-flip field decreases and this is reflected well in the MR data of EuIrAl 4 Si 2 at 2, 5 and 10 K, respectively. Just above T N1 at 15 K the MR of EuIrAl 4 Si 2 in the paramagnetic state is relatively small and negative in almost the entire range up to 8 T. The negative, albeit small, MR at 15 K may be due to the combined partial quenching of the spin fluctuations above T N1 and positive MR due to background conduction electrons.

The MR for H ‖ [001] and J ‖ [100] at selected values of temperature in fields up to 8 T is shown in Fig. 4(f). The prominent feature is the sharp drop and hysteresis in MR at fields where the metamagnetic jump is seen in the magnetization of the Ir compound. The hysteresis is reduced in width and occurs at lower fields with the increase in temperature, in excellent correspondence with the magnetization data for H ‖ [001] (see, Fig. 2(b) inset). Thus there is a close correlation between the magnetization and the MR data which clearly reveal the first order nature of the metamagnetic transition. Qualitatively the behavior of MR in the two compounds is similar. At low fields, below T N , it is positive and increases with field. Above the metamagnetic transition the MR increases with field, its magnitude and sign depending upon the temperature and field. In the paramagnetic state at 15 K in EuRhAl 4 Si 2 and at 20 K in EuIrAl 4 Si 2 , the MR is negative at all fields.

151 Eu Mössbauer spectra. Additional information on the two compounds was obtained by employing 151 Eu Mössbauer spectroscopy. Spectra were recorded between 4.2 and 15 K. We show here only the data for T = Ir (the spectra for the Rh-analogue are shown in the Supplementary information). We recall that the specific heat data show a transition near 14.7 K followed by another one at 13.2 K. At 4.2 K, the spectrum is a standard Eu 2+ hyperfine pattern with a single hyperfine field of 30.1(1) T, meaning that the magnetic structure is equal moment. The spectra start changing shape near 13 K, and Fig. 5 shows their rapid evolution in the small temperature interval 13-15 K. At 13, 14 and 14.5 K, the spectra are purely incommensurate modulation (ICM) patterns, and at 15 K the spectrum is a superposition of an ICM spectrum (red line) and of a single line paramagnetic spectrum (green line). Above 15 K, the spectrum is a single line characteristic of the paramagnetic phase. Although the transitions we observe by Mössbauer spectroscopy are slightly shifted with respect to those given by the specific heat peaks, the spectra illustrate well the way EuIrAl 4 Si 2 goes, on cooling, from the paramagnetic phase to the ICM phase, then to the commensurate lock-in phase, the first transition being first order as witnessed by the coexistence of the two phases in the spectra at 15 K. The spectra in EuRhAl 4 Si 2 are very similar, the transitions being shifted to somewhat lower temperature, in agreement with the specific heat data indicating a first transition close to 11.7 K, followed by another at 10.4 K. This cascade of transitions is a rather common feature in intermetallics with divalent Eu or trivalent Gd ions [START_REF] Blanco | Specific heat in some gadolinium compounds. II. Theoretical model[END_REF] . In particular, it is present in EuPdSb [START_REF] Bonville | Incommensurate modulated magnetic structure in orthorhombic EuPdSb[END_REF] and has been discovered more recently in EuPtSi 3 [START_REF] Kumar | Magnetic properties of EuPtSi 3 single crystals[END_REF] and EuNiGe 3

2

. It seems to be concomitant with the occurence of a sizeable anisotropy of the magnetisation, as in the present materials.

Discussion

Two noticeable features in the magnetization along the c-axis in both compounds are the existence of a jump at low field followed by a plateau region with magnetization 1/3 M 0 (M 0 = 7 μ B /f.u. at 2 K). In some rare earth antiferromagnets, plateaus in the magnetization along the easy-axis have been reported, but at a finite field value due to the strong crystal field induced Ising character of the considered rare earth [START_REF] Gignoux | Magnetic properties of intermetallic compounds[END_REF] . Some examples are DyCo 2 Si 2 (T N = 24 K) [START_REF] Iwata | Metamagnetism in DyCo 2 Si 2[END_REF][START_REF] Shigeoka | Magnetic Phase Diagram and Field Induced Magnetic Structure in DyCo 2 Si 2[END_REF] , orthorhombic DyCu 2 (T N = 31 K) [START_REF] Iwata | Magnetization process of single crystal DyCu 2[END_REF] and CeZn 2 (T N = 7 K) [START_REF] Gignoux | Magnetic phase diagrams and metamagnetic processes in Ising systems: The case of CeZn 2[END_REF] . The intermediate phases at the plateaus correspond to successive flipping of blocks of magnetic moments which change the magnetic periodicity and the propagation vector. In hexagonal PrRh 3 B 2 , which shows the intermediate 1/3M 0 state, it was proposed that the ferromagnetic linear chains along [0001] are coupled antiferromagnetically, which leads to an Ising-spin structure in a triangular lattice [START_REF] Yamada | Magnetic and Fermi surface properties in PrRh 3 B 2[END_REF] .

To our knowledge, the only very low field 1/3 M 0 plateau has been observed at very low temperature for a field along [111] in spin-ice Dy 2 Ti 2 O 7 , which does not undergo magnetic ordering. In this case, the plateau is due to the peculiar pyrochlore lattice structure with extreme Ising anisotropy at each of the 4 Dy sites having threefold symmetry along one of the 111 axes.

In EuRhAl 4 Si 2 and EuIrAl 4 Si 2 , the presence of such an unusual low field plateau is probably due to a combination of a low (but sizeable) crystalline anisotropy and of a peculiar ferrimagnetic-like magnetic structure, which we believe consists in a tripled magnetic unit cell along c with "up-up-down" moment arrangement. The propagation vetor would then be k = (0 0 1/3). Neutron diffraction experiments on a single crystal of EuRhAl 4 Si 2 are planned in order to check the existence of such a structure.

We present in the following a model which sets our assumption of an "up-up-down" sequence, with moments along c, on a firmer ground. The model considers a magnetic super-cell with 6 sub-lattices consisting each of a ferromagnetic plane perpendicular to [001]. The planes are bi-dimensional simple cubic lattices with parameter a, and the distance between planes is c. The Eu ions are taken to interact via the short range RKKY exchange and the infinite range dipolar interaction, and they are submitted to a crystalline anisotropy of the form:

= H DS an z 2
, where Oz is the c-axis, appropriate for axial symmetry around c. As underlined in Refs 2,13, inclusion of the dipolar interaction is essential to reproduce the magnetic behaviour in a realistic way in compounds with S-state Eu 2+ and Gd 3+ ions because it competes with the crystalline anisotropy for establishing the moment direction and the spin-flop fields. We detail in the following these three interactions: i) we consider that each Eu ion interacts with its neighbours through 3 exchange integrals: in-plane J 0 (4 nearest neighbour ions at distance a), and interplane (along c) J 1 (2 first neighbour ions at distance c) and J 2 (2 second neighbour ions at distance 2c). We treat the exchange in the molecular field approximation, i.e. the exchange field H ex i on a given spin S i writes:

+ - + - g J J J H 1 [4 S S S S S ] 2 ex i B i i i i i 0 1 1 1 2 2 2
where the indexes must be taken mod. 6. Due to the oscillatory character of the RKKY interaction, the integrals can be positive (ferromagnetic) or negative (antiferromagnetic).

ii) the infinite range dipolar field on each ion contains two contributions: the field originating from the ion own sub-lattice and that coming from each of the 5 other sub-lattices. The field acting on an ion in sublattice i can be written:

= ∑ , H Q m dip i k k i k
where Q k i is a matrix containing lattice sums evaluated using Ewald-type summations [START_REF] Wang | Estimate of the cutoff errors in the Ewald summation for dipolar systems[END_REF] , k runs over all the sublattices and m k is the moment in the kth sublattice (in units of μ B ). In units of the Lorentz field = in the same units. For EuRhAl 4 Si 2 , a = 0.418 nm and c = 0.828 nm so that H L = 44.76 G. Since m 0 = 7 μ B , the dipolar fields are of the order of 0.3 T.

iii) the crystalline anisotropy term DS z 2 would create, in the absence of dipolar field, an easy axis along [001] (c-axis) for D < 0. In the presence of the dipolar field, we find that D must be taken rather large (negative) in order to have moments along c, otherwise the dipolar field forces the moment to lie in the (ab) plane.

In order to obtain the magnetisation and the susceptibility, we perform for each sub-lattice an exact diagonalisation of the hamiltonian:  ( )

μ = + . + + , ( ) DS g S H H H 4 i z B dip i ex i 2
where H is the applied magnetic field, and obtain S i which is treated self-consistently within the 6 sublattices until convergence is reached. A scan in the 4 parameter space {J 0 , J 1 , J 2 , D} was performed, taking the nearest neighbour in-plane exchange integral J 0 positive (ferromagnetic), ensuring ferromagnetic planes normal to [001], and D large negative, ensuring the c-axis as easy axis. We could find for both compounds a point around which the zero field ground state is the "up-up-down" structure, and for which the magnetisation curves are in rather good agreement with experiment as to the low field jump, the spin-flop and spin-flip fields. However, the T N values are somewhat lower (by 2-3 K) than the actual values. For EuRhAl 4 Si 2 , these parameters are: J 0 = 0.18 K, J 1 = -0.02 K, J 2 = -0.17 K and D = -1.3 K. The corresponding m(H) and χ(T) curves are shown in Fig. 6(a), where it can be seen that .  T 7 5K

N , instead of the experimental value 10-11 K. Similarly, the parameter set: J 0 = 0.22 K, J 1 = -0.025 K, J 2 = -0.29 K and D = -1.5 K reproduces well the data in EuIrAl 4 Si 2 (see Fig. 6(b)), but with T N = 10 K instead of 13-14 K. For this latter compound, J 0 and J 2 are larger than for EuRhAl 4 Si 2 , which accounts for the higher Néel temperature and characteristic field values. It may be noted that J 1 is close to zero, which could correspond to a node of the RKKY oscillating function, and that the crystalline anisotropy coefficient D is unusually large in these two materials. Indeed, although this quantity is seldom measured, its expected order of magnitude for Eu 2+ or Gd 3+ is 0.1 K 15 . We presently have no explanation for this enhanced value.

Conclusion

We have performed a series of measurements to probe in detail the magnetic properties of single crystals of two tetragonal quaternaries, EuRhAl 4 Si 2 and EuIrAl 4 Si 2 . Both present a first order transition from a paramagnetic to an incommensurate antiferromagnetic phase at 11.7 and 14.7 K respectively, followed by a lock-in transition to a commensurate antiferromagnetic phase at 10.4 and 13.2 K. Such a transition cascade is rather common among Eu intermetallics. By contrast, the magnetic behavior deep in the antiferromagnetic phase is rather unexpected for isotropic divalent Eu: the response is antiferromagnetic for a field applied in the ab-plane and ferromagnetic for a field applied along c. In this latter case, the magnetization shows a low field jump followed by an intermediate plateau region attaining a value of 1/3 M 0 . Beyond the plateau region there is a sharp hysteretic spin-flip transition. Our in-field transport data correlate well with this behaviour of the magnetisation. Surprisingly, the magnetization shows a substantial hysteresis around the origin which is not expected in divalent Eu compounds. The magnetization data hint to the presence of an unusual "up-up-down" zero field magnetic structure which is shown to be obtainable by the mean field model presented here, but which must be confirmed by neutron diffraction experiments.

Methods

The details of the single crystal growth and structure have been reported in Ref. 1. In the present work single crystals of LaTAl 4 Si 2 (T = Rh and Ir), isostructural to Eu analogs, were also grown following the same protocol as adopted for the Eu compounds. The flux-grown crystals were oriented by means of the Laue diffraction technique in back-reflection mode and cut appropriately using a spark erosion cutting machine. The magnetization as a function of temperature and magnetic field was measured in Quantum Design SQUID and VSM magnetometers between 1.8 and 300 K. Heat capacity and electrical resistivity in zero and applied fields were measured in a Quantum Design PPMS. 151 Eu Mössbauer absorption spectra were taken at selected temperatures using a commercial 151 SmF 3 γ-ray source, mounted on a constant acceleration electromagnetic drive, in a standard liquid He cryostat. This work is licensed under a Creative Commons Attribution 4.0 International License. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in the credit line; if the material is not included under the Creative Commons license, users will need to obtain permission from the license holder to reproduce the material. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/

  Magnetic anisotropy, unusual hysteresis and putative "up-updown" magnetic structure in EuTAl 4 Si 2 (T = Rh and Ir) Arvind Maurya 1 , A. Thamizhavel 1 , S. K. Dhar 1 & P. Bonville 2

Figures. 1

 1 (a,b) show the inverse susceptibility data between 300 and 1.8 K respectively for EuRhAl 4 Si 2 and EuIrAl 4 Si 2 in a field of 0.1 T applied along the [001], [100] and [110] directions. The molar susceptibility follows the Curie-Weiss law:

Figure 1 .

 1 Figure 1. Inverse magnetic susceptibility of (a) EuRhAl 4 Si 2 and (b) EuIrAl 4 Si 2 along [001], [100] and [110] crystallographic directions. The insets show the anisotropy in the low temperature susceptibility. (c-d) The magnetic susceptibility in a field of 0.005 T: zero field cooled (ZFC), field cooled cooling (FCC) and field cooled heating (FCH) data for H ‖ [100] and H ‖ [001].

  Figures 4(a,b) depict the variation of electrical resistivity ρ(T) with temperature, which reveals a normal metallic nature for the two compounds. The data have been taken with the current density along [100] and [001], respectively. The resistivity is anisotropic which indicates an anisotropic Fermi surface of these two compounds with tetragonal symmetry.

Figure 2 .

 2 Figure 2. M vs H for (a) EuRhAl 4 Si 2 and (b) EuIrAl 4 Si 2 at 1.8 K for H in the ab-plane and along [001] and at selected temperatures for H along [001](insets) (c). Low field hysteresis showing the coercive field in EuRhAl 4 Si 2 (c) and EuIrAl 4 Si 2 (d).

Figure 3 .

 3 Figure 3. Heat capacity C p and entropy S mag as a function of temperature in EuTAl 4 Si 2 and nonmagnetic LaTAl 4 Si 2 for (a) T = Ir and (b) Rh. The horizontal line depicts the high temperature limit (R In 8) of the magnetic entropy for an S = 7/2 spin system.

Figure 4 .

 4 Figure 4. Electrical resistivity as a function of temperature when the current is parallel to major crystallographic directions in (a) EuRhAl 4 Si 2 and (b) EuIrAl 4 Si 2 . ρ(T) data as a function of temperature at selected magnetic fields applied along (c) [010] and (d) [001], and field evolution of MR for (e) H ‖ [010] and (f) H ‖ [001] in EuIrAl 4 Si 2 with current being parallel to [100].

Figure 5 .

 5 Figure 5. 151 Eu Mössbauer spectra in EuIrAl 4 Si 2 between 4.2 K and 15 K (see text) . The modulations at 13 K, 14 K and 14.5 K along the propagation vector k are shown next to the corresponding spectrum.

3 ,

 3 with V = 6a 2 c being the volume of the chosen magnetic cell, the own ion dipolar matrix writes: matrices are equal with only one non-zero matrix element:

Figure 6 .

 6 Figure 6. Calculated isothermal magnetization curves at 1.8 K for H || [001] and in the (ab) plane, and, in the insets, calculated susceptibility for H || [100] : (a) with J 0 = 0.18 K, J 1 = -0.02 K, J 2 = -0.17 K and D = -1.3 K, to be compared with data in EuRhAl 4 Si 2 of Figs. 1(a) and 2(a) and (b) with J 0 = 0.22 K, J 1 = -0.025 K, J 2 = -0.29 K and D = -1.5 K, to be compared with data in EuIrAl 4 Si 2 of Figs. 1(b) and 2(b).

  How to cite this article:Maurya, A. et al. Magnetic anisotropy, unusual hysteresis and putative "upup-down" magnetic structure in EuTAl 4 Si 2 (T = Rh and Ir). Sci. Rep. 5, 12021; doi: 10.1038/srep12021 (2015).

Table

EuRhAl 4 Si 2 EuIrAl 4 Si 2 μ eff

	H ‖ [100]	7.64	8.0	7.78	3.0
	H ‖ [110]	7.84	6.9	8.06	2.3
	H ‖ [001]	7.83	10.4	7.90	9.2

4

Si 2 , (μ B /f.u.) θ p (K) μ eff (μ B /f.u.) θ p (K)

Table 1 . Effective moment and paramagnetic Curie temperature values in EuRhAl 4 Si 2 and EuIrAl 4 Si 2

 1

T N1 T N2 Spin flip field at 1.8 K (T)

					Hysteresis at 1.8 K
						H ‖ c
						Plateau
			H ‖ ab-plane	H ‖ c	H c (T)	width (T)
	EuRhAl 4 Si 2	11.7 10.4	5.5	1.8	0.18	1.2
	EuIrAl 4 Si 2	14.7 13.2	9.2	3.9	0.25	2.5

Table 2 . Magnetic characteristics of EuRhAl 4 Si 2 and EuIrAl 4 Si 2 .

 2

Scientific RepoRts | 5:12021 | DOi: 10.1038/srep12021 www.nature.com/scientificreports

Scientific RepoRts | 5:12021 | DOi: 10.1038/srep12021

Author Contributions

S.K.D. and A.T. proposed the single crystal growth and the study of their anisotropic magnetic properties. A.M. has grown the crystals and did susceptibility, heat capacity and electrical transport measurements, and analyzed the data. Mössbauer spectra were recorded and analyzed by P.B. The main manuscript was written by S.K.D., P.B. and A.M. All authors have reviewed the manuscript.

Additional Information

Supplementary information accompanies this paper at http://www.nature.com/srep Competing financial interests: The authors declare no competing financial interests.


Osvaldo L Asato
email: osvaldo.asato@usp.br

Guilherme M Dobrianskyj

Fabrício Junqueira

Diolino J Santos Filho
email: diolinos@usp.br

Paulo E Miyagi
email: pemiyagi@usp.br

Process Control System Considering the Machines Functional Flexibilities

Keywords: Multifunctional Machine Tools, Functional flexibility, Petri net, Ant Colony Optimization

The evolution of tool-machines shows a trend for integration of different machining functionalities such as, milling, turning, and drilling, on a same machine. These machines are result of the mechatronics development that allows the aggregation of different functionalities into a unique machine (named multifunctional machine tool -MMT). However from the point of view of the production process there is a limitation of the conventional process control solution, i.e., when the systems are composed by MMTs, their functional flexibilities are considered only during scheduling and encapsulated at this level, and during the process execution, MMTs are viewed as traditional ones. Then, the multi-functionalities are not explored for resources allocation. Thus, a new process control approach is proposed here to dynamically allocate the resources exploring the functional flexibilities of MMTs.

Introduction

The technological developments in mechatronics systems have allowed the aggregation of several machining functionalities (such as milling, turning and drilling) at a same machine tool [START_REF] Moriwaki | Multi-functional machine tool[END_REF]. These multifunctional machine tools (MMTs) have then a large functional flexibility than conventional machine tools. These MMTs potentially meet the need for productivity and efficiency of production processes, especially in the production of customized products [START_REF] Aronson | Machine-tools basics part1[END_REF], [START_REF] Shin | Self-evolution framework of manufacturing systems based on fractal organization[END_REF]. However, currently the process control solutions for productive systems (PSs) do not consider and do not explore any nondeterministic behavior that the MMTs flexibility introduces in these systems. Conventionally, MMTs are used as a type of machining center [START_REF] Wiendalhal | Changeable manufacturing classification, design and operation[END_REF], where the focus is on the centralization of the functionalities in the same machine (without work piece manipulation between different machines). In this context, a deterministic (static) allocation of resources (machine tools, machining centers, MMTs) is adopted according to the production plan [START_REF] He | Research on job shop scheduling aiming at instant customerisation[END_REF], [START_REF] Li | A review on integrated process planning and scheduling[END_REF], and in the case of intrinsic indeterminism, such as delays or early finish of production activities, the potential flexibility of MMTs are not used to maintain or improve the PS's global performance.

Therefore, to explore the functional flexibility of MMTs, a new process control solution must be considered. The focus must be on both functional and operational flexibilities of the MMTs when the resources are designated in real time according to the functionality demand from the production process. This approach changes the paradigm of the resource allocation control that conventionally is based on the predefined functionality of the available resources. The proposed control solution for PS considers the necessary functionalities for the execution of each activity of production process, and establishes a dynamic allocation of resources. The adopted modeling tool considers that production processes and its control can be treated as a class of discrete event system [START_REF] Cangea | Designing a control structure for discrete event systems described by Petri nets[END_REF], [START_REF] Matsusaki | Modeling of distributed collaborative control systems of production systems[END_REF], [START_REF] Santos Filho | Task Control of Intelligent Transportation Vehicles In Manufacturing Systems[END_REF], [START_REF] Nakamoto | Resources allocation control in flexible manufacturing systems using the deadlock avoidance method[END_REF]. This paper is organized as follows: Contribution to value creation in section 2; Fundamental concepts in section 3; Process Control with functional flexibility in section 4; Case study in section 5 and the Conclusion in section 6.

Contribution to Value Creation

Although the technological evolution has developed resources with multiple functionalities as MMTs, it wasn't developed a control system that, at the moment of the allocation of the resource the shop floor, which of the functionalities of this machine is the most appropriate to use in relation the global production process.

In current Resource Allocation Control Systems (RACS) for PS, MMTs are used as conventional machine, because the control system works with static allocation, in other words, the functionalities are attributed for a local processing in such a way that RACS doesn't have access to this type of information. In this context, the current concept of RACS doesn't explore the potentiality of MMTs, once the same was not projected to allocate functionalities. It is observed that the solutions of RACS use graphs of allocation of resources [START_REF] Nakamoto | Resources allocation control in flexible manufacturing systems using the deadlock avoidance method[END_REF], but there are not proposals in the sense of building, in equivalent way, graph of allocation of functionalities. In fact to explore the potentiality of MMTs is necessary to modify the structure of RACS, due to emergence of one more variable control called functionality. In the conventional control systems there were a one to one relationship between resource and functionality, sparing additional control strategies. However, if this reality was modified in which a resource possess a group of associated functionality, is necessary to establish specifics control rules to treat this new variable. Taking in account a new paradigm is proposed in terms of modular structure for a in that the control of allocation of functionalities to the resources in real time it is modeled to improve the efficiency of the global system [START_REF] Asato | Control of Productive Systems with Functional Flexibility Level[END_REF]. The term efficiency should be understood as smaller cost for production: (i) it can decrease the time of execution; (ii) it can avoid that the system enters in deadlock for not being capable of reconfiguring the plan of the allocation of resources.

The obtained result will be: (i) to perfect the rational use of the resources in productive systems; (ii) to increase the competitiveness among the companies; (iii) to generate benefits for the society that can have access to products of lower cost.

Fundamental Concepts

The production process as well as its control procedure can be treated as a processes with discrete states that are abruptly changed when an instantaneous event occurs, for example machine that no working. This process defines a class of discrete event system [START_REF] Cangea | Designing a control structure for discrete event systems described by Petri nets[END_REF], [START_REF] Matsusaki | Modeling of distributed collaborative control systems of production systems[END_REF], [START_REF] Santos Filho | Task Control of Intelligent Transportation Vehicles In Manufacturing Systems[END_REF], [START_REF] Nakamoto | Resources allocation control in flexible manufacturing systems using the deadlock avoidance method[END_REF] and therefore, modeling tools such as Petri net [START_REF] Cangea | Designing a control structure for discrete event systems described by Petri nets[END_REF] can be used not only to describe the structure and behavior of the processes, but also as a specification of a desired task execution.

Based on previous works related with PSs, the Enhanced Mark Flow Graph (E-MFG) [START_REF] Santos Filho | Task Control of Intelligent Transportation Vehicles In Manufacturing Systems[END_REF] is adopted to model the production processes, and to specify its control solution. The E-MFG is derived from Petri net considering individualized marks through the association of attributes. With E-MFG is possible to control the information flow in the model, i.e., modifying the attributes of the marks and, based on this information, designate resource using an element denominated control-box. In this context, the control-box works with ant colony algorithm to designate the machine in the shop floor.

Enhanced -Mark Flow Graph (E-MFG)

The E-MFG is effective to model control strategies, especially when associated to the concept of "functional blocks" for modular and distributed control as established in the standard IEC 61499. From the practical point of view of SPs, marks in E-MFG are associated to a vector of attributes that guarantees their individuality. These attributes can be associated to the information about products, processes, transportation resources, transformation resources etc.. The manipulation of the mark's attributes indicates the alteration of some information that is interpreted as an updating of the local states of the system (state associated with a box that represents an activity, operation or task of the production process). The control-box accomplishes this function based on ant colony algorithm that is applied for the verification and update the attributes previously specified, and also according to the firing rule of the transitions, as exemplified in Fig. 1. In this example, a mark with the attributes <a1,a2,a3,a4> where: a1=functionality, a2= type of work-piece, a3= origin (from where the work-piece is coming from), a4= allocation (where the piece should go).In PS, there is an intense interaction among the entities that compose the system [START_REF] Matsusaki | Modeling of distributed collaborative control systems of production systems[END_REF] and E-MFG has also has a specific graph structures called transmission and reception interface [START_REF] Nakamoto | Resources allocation control in flexible manufacturing systems using the deadlock avoidance method[END_REF].

The transmission interface sends an asynchronous message when the respective box is marked. The transmission interface corresponds to the "method call" of an "object" based on the object-oriented paradigm. This structure models the message emission for a system and it is composed by a send transition and a send arc (Fig. 2a). The reception interface enables or inhibits the transition firing connected to it, in order to the received information. It models the capture and the processing of messages coming from other systems. This structure is composed of a receive box and a receive arc (Fig. 2b).

Algorithm Ant Colony Optimization

The algorithm "ant colony optimization" (ACO) was developed by Dorigo [START_REF] Dorigo | The ant colony optimization metaheuristic[END_REF] and in his work; this algorithm was used to analyze the problem of resource designation. In this context the ACO verifies which resources have the functionality required and among them the resource that has a functionality with lower operating cost is selected.

The state transition rule is as follows:

i. an ant positioned on node r chooses the resource s to move to by applying the rule given by Eq. ( 1).

Arg max { [τ(r,j)] α . [η(r,j)] β }; jϵ Cn(r), se q ≤ qo S = (Eq.1) S, otherwise.

Where: S is the next resource, r = previous resource, n = total number of resources, Cn(r) = group of the resources that possess the following functionality, S = defined resource following the distribution of probability p(r,s).

ii. P(r,s) probability for an ant to travel from node r to node s, the rule given by Eq.( 2).;

[τ(r,s)] α . [η(r,s)] β if s ϵ Cn(r) Ʃ iϵcn(r) [τ(r,i)] α . [η(r,i)] β (Eq.

2) 0 otherwise

Where: n = total number of resources, Cn(r) = group of the resources that possess the following functionality. Pheromone level based on local updating rule: τ(r,s)←(1-ρ).τ(r,s)

iv.

Pheromone level based on global updating rule: τ(r,s)←(1-α).τ(r,s)+α.∆τ(r,s)

Process Control with Functional Flexibility

The approach adopted is based on the semantic division of the information considered for the production process control: (a) the information related to functionalities requested to execute an activity, and (b) the information related to allocation of resources in accordance with the functionality requested. This approach allows resources to be designated in real time according to the functionality demand.

RACS Structure

In this work is proposed an evolution [START_REF] Asato | Control of Productive Systems with Functional Flexibility Level[END_REF] according to Fig. 3a where exists: (i) a specific module for production sequence control (PsC); (ii) a module for functionality designation (FD) according to the set of activities to be executed; and (iii) a module for resources designation based on the functionalities (RD-F). Each module accomplishes specific control functions and, through the cooperation among them, production activities (in accordance with a global production plan) are executed, as described to proceed:

• PsC represents the sequence of operations to be executed for the production of one product. This control module needs to request resources, and then receives information about the designation of resources to be allocated for the execution of the set of the operations in subject. • The FD controls module identifies the functionality associated with the requested resources from PsC, according to the original production plan. This module specifies the functionality corresponding to the resources requested by the PsC. • The RD-F controls module designates the resources, i.e., RD-F determines the association between functionality and the best resource to effectively execute each process operation.

RD-F Algorithm

The resource designation is made based on ACO that is processed whenever there is a transition of states in the FD control module: the RD-F module considers the actual state of each machine including MMTs, and also its operational and functional flexibility. These modules interact in the following way: (i) PsC monitors the evolution of the production processes and when a new activity must be executed, it sends a request to FD about the needed resources according to the original production plan; (ii) based on this request, FD identifies the respective functionality and sent this information to RD-F that accomplishes the designation of the resource that can execute the requested functionality, according to ACO: the solution corresponds to the vector of Resources that possess the largest pheromone. This solution is sent back for PsC; (iii) Once this information is received, PsC allocates the specified resource, and control the execution of the activity production.

Synthesis of the Control Modules

The control solution, i.e., the procedures to be executed in each module of RACS can be specified through E-MFG models and ACO algorithm. Fig. 4 illustrates the procedure to develop the models used for decision making. The procedure is composed of the following steps:

Step 1 -FSG Synthesis. Based on each production process model, a general structure of "functionalities sequence graph" (FSG) is developed in E-MFG. Each resource is mapped according to the functionality that is necessary to execute the corresponding activity (see Fig. 4).

Step 2 -Refinement of FSG. For each activity of the production process, FSG identifies a specific functionality. The information about the resource can be encapsulated through the concept of individual marks of the E-MFG. Therefore, the vector of mark attributes is defined as: <a1,a2> where, a1=Fi (functionality Fi to execute an activity of the production process) and a2=Ri (resource Ri according to the production process model previously specified). In this context, the graph structure between transitions Ai and Ai+1 of FSG is refined trough a graph with a box (input buffer), a control-box (representing the functionality) and a box (exit buffer) (Fig. 5). The rules associated with the control-box will be defined in step 4. Step 3 -FD Synthesis. From refined FSG of each production process, the elements with the same semantics are identified and merged, ie, control-boxes related to the same functionality are fused. The result E-MFG global model of all refined FSG is the control model for FD module. The problem of conflict when there is dispute of resources for a same functionality is solved automatically looking for the solution of smaller global cost according to ACO.

Step 4 -RD-F Rules Synthesis. In this step, the ACO is applied to generate the best designation of resources (lower operation cost) whenever there is a state transition in FD. The cost table associated to ACO is based on the information defined previously in the scheduling phase. Therefore, the rule for resources designation is executed based on ACO and the solution of the algorithm is used to update the mark's attributes of the E-MFG model of FD.

Step 5 -PsC Synthesis. It is derived from control models for FSG as is it follows:

• It is added a complementary box related to each box Fi of the FSG models.

• Each box Fi between transitions Ai and Ai+1 of FSG is substituted by a graph with: (i) a box to model the request of a resource, (ii) a box to represent the allocation of a resource and (iii) a box to indicate that the resource is available.

• The interaction between PsC, FD and RD-F modules are established through communicator arcs (send and receive arcs) Fig. 3b.

Case Study

Consider a PS that has a production process with:

• Functionalities F1, F2, F3.

• Machines R1, R2, R3, R4, R5, R6. Each resource has the functionalities according to Table 1. Hypothesis: • All machines are available at time t=0. The machines do not break.

• The transportation time between machines is not considered

• The setup time is not considered for all machines.

The production process has the following the sequence: F1→ F2→ F3.

Taking in account the worst case is when the functionalities are being requested simultaneously. The Table 2 represents the costs associated, and the first column (resources R1, R2, R3, R4, R5, R6) refers to the machine of origin (start) and the first line refers to the machine of destiny (next machine allocated). When beginning the process in one of the machines, an initial cost exists (start of the process) according to Table 2.

Table 2.. The initial cost attributed to the resources.

In same way for the last machine Table 3. When processing the ACO algorithm was considered: • The following parameters: α = 1; β = 2; ρ = 0.9; ϕ = 0.2; qo = 0.5.

• Initial value of the pheromone = 0.3. The solution is: R2 -R6 -R5.

Conclusion

This approach uses algorithm ACO, for each scenario were performed 50 iterations on Core™ 2 Quad CPU Q9550 2.84Ghz 4GB of RAM and Windows 7 system. In this condition the algorithm takes 16 milliseconds to meet the result (R2, R6, R5) with a cost 9 units. Therefore, the proposed approach for the production process control system using functionality allocation makes possible the dynamic allocation of resources in agreement with the requested functionality of the process. This aspect is fundamental to explore the potential of production systems (PSs) with multifunctional machine tools (MMTs). The process control system proposed in this work can dynamically make decisions in accordance with the monitored signals of the PS exploring the flexibility of the available resources when activity delays of early finishing occurs. These are the first result using ACO. The following next researches are about the use of this approach to avoid deadlock in complex PS.

Fig. 1 . 2 .

 12 Fig. 1. E-MFG with transitions firing Fig. 2. Example of communicators.

( b )

 b Interface of reception (a) Interface of send P(r,s)= Process Control System Considering the Machines Flexibilities 137 iii.

Fig. 3a . 4 .

 3a4 Fig.3a. Modules of the control Fig.3b. communication between Fig.4. Procedure to synthesis system modules the models.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Control model for FSG. Fig. 6. Refinement of control model for FSG.

Table 1 .

 1 Description of the functionalities of each resource of the PS.

	Resource functionality Resource functionality	Resource functionality
	R1	F1	R3	F2	R5	F2 and F3
	R2	F1	R4	F2	R6	F2 and F3

Table 3 .

 3 The end cost attributed to the resources.

Table 4 .

 4 Table Cost associated to the resources.

	From	To	R1	R2	R3	R4	R5	R6
	R1		1	3	2	8	4	5
	R2		3	1	5	13	8	2
	R3		3	4	1	7	4	11
	R4		10	8	7	1	5	3
	R5		7	8	1	5	1	4
	R6		1	3	2	7	3	1

•

  Number of iterations: 50 (with 10 ants). The result is presented in the Table5.

	R1=2.5	R2=2	R3=3 R4=10	R5=7 R6=6
	R1=10	R2= 9	R3=3	R4= 4	R5= 2	R6= 5

Table 5 .

 5 Obtained solution using algorithm based on ACO.

	Ant	F1	F2	F3	Cost	Ant	F1	F2	F3	Cost
	1	R2	R6	R6	10	27	R2	R6	R6	10
	2	R1	R3	R5	10,5	28	R1	R3	R5	10,5
	3	R2	R3	R6	23	29	R2	R6	R6	10
	4	R1	R5	R5	9,5	30	R1	R3	R6	15,5
	5	R2	R6	R5	9	31	R2	R6	R5	9
	6	R1	R4	R6	18,5	32	R2	R6	R5	9
	7	R2	R4	R5	22	33	R2	R6	R5	9
	8	R1	R6	R6	13,5	34	R2	R6	R5	9
	9	R2	R5	R6	19	35	R2	R3	R5	13
	10	R1	R6	R5	12,5	36	R2	R6	R6	10
	11	R2	R6	R5	9	37	R2	R4	R5	22
	12	R2	R6	R6	10	38	R2	R3	R5	13
	13	R2	R3	R5	13	39	R2	R5	R6	19
	14	R2	R5	R5	13	40	R1	R6	R5	12,5
	15	R2	R4	R5	22	41	R2	R6	R5	9
	16	R2	R6	R5	9	42	R2	R6	R5	9
	17	R1	R4	R6	18,5	43	R2	R6	R5	9
	18	R2	R3	R6	23	44	R2	R6	R5	9
	19	R1	R5	R5	9,5	45	R2	R6	R5	9
	20	R2	R6	R6	10	46	R2	R5	R5	13
	21	R2	R6	R5	9	47	R2	R4	R6	23
	22	R2	R6	R5	9	48	R2	R4	R6	23
	23	R2	R5	R6	19	49	R2	R4	R5	22
	24	R2	R6	R5	9	50	R2	R3	R6	23
	25	R2	R4	R6	23
	26	R2	R5	R5	13

Process Control System Considering the Machines Flexibilities 135

Acknowledgments. The authors thank the financial support of the Brazilian government agencies, CNPq, FAPESP and to the MEC/CAPES/PET program.