D Truzzolillo
email: domenico.truzzolillo@univ-montp2.fr

V Roger

C Dupas

S Mora

L Cipelletti

Bulk and interfacial stresses in suspensions of soft and hard colloids

Keywords: numbers: 82.70.Dd, 68.05.-n, 83.80.Hj

We explore the influence of particle softness and internal structure on both the bulk and interfacial rheological properties of colloidal suspensions. We probe bulk stresses by conventional rheology, by measuring the flow curves, shear stress vs strain rate, for suspensions of soft, deformable microgel particles and suspensions of near hard-sphere-like silica particles. A similar behavior is seen for both kind of particles in suspensions at concentrations up to the random close packing volume fraction, in agreement with recent theoretical predictions for sub-micron colloids. Transient interfacial stresses are measured by analyzing the patterns formed by the interface between the suspensions and their own solvent, due to a generalized Saffman-Taylor hydrodynamic instability. At odd with the bulk behavior, we find that microgels and hard particle suspensions exhibit vastly different interfacial stress properties. We propose that this surprising behavior results mainly from the difference in particle internal structure (polymeric network for microgels vs compact solid for the silica particles), rather than softness alone.

Introduction

Colloidal suspensions comprise solid particles in the size range 1 nm -10 µm suspended in a background fluid. They are ubiquitous in every day life and industrial applications, and they are intensively studied as model systems in condensed matter, e.g. for tackling problems such as the glass transition or crystallization. Colloidal suspensions, even in the concentrated regime or in the presence of strong interparticle interactions, are soft systems, in that they are significantly deformed or even driven from a solid-like to a fluid-like state by relatively modest external forces. Indeed, how colloidal systems may respond to an external drive is a key property in many applications [START_REF] Larson | The structure and rheology of complex fluids[END_REF], besides being a rich and fascinating problem per se. Accordingly, there has been a great interest in the rheological properties of suspensions, both in the linear regime corresponding to small drives and in the non-linear regime [START_REF] Koumakis | [END_REF]3]. The availability of particles in a large variety of shapes [4,5,6,7,8,9,10], with different degrees of softness [START_REF] Koumakis | [END_REF]11,12,13], and with highly tailored interactions (including selective interactions mediated by DNA [14] or anisotropic interactions as in patchy colloids [9]) has further spurred studies on the relationship between the rheological properties of a suspension and the structure and interactions of its constituents. In particular, the role of particle softness has been investigated in several recent works, in part thanks to the availability of particles with a controlled degree of softness, e.g. solid particles covered by a polymer layer [START_REF] Koumakis | [END_REF]15] or microgel particles formed by cross-linked polymers [11,16]. There is some consensus that the rheological properties of small, submicron hard particles are governed by thermal stresses of entropic origin, which sharply increase well before attaining the random close packing volume fraction, while the flow properties of jammed packings of large, athermal and deformable particles are governed by the elastic energy associated with particle deformation. The behavior of particles at intermediate size and softness scales, such as sub-micron PNiPAM microgels, however, is much more controversial and the nature of the fluid-to-solid transition (glass transition vs jamming) in these systems is highly debated [3,17,18,19]. Soft matter systems may also exhibit intriguing behaviors ruled by interfacial stresses, as opposed to the bulk stresses discussed above. For example, recent work has shown that the addition of colloidal particles at the interface between two (immiscible) fluids may modify profoundly the mechanical properties of the interface, by imparting it mechanical rigidity as in "bijels" [START_REF] Stratford | Colloidal jamming at interfaces: A route to fluidbicontinuous gels[END_REF]21] or in non-spherical "armoured" bubbles [22,23], or by protecting bubbles and drops from coalescence [24]. Recently, we have explored the role of interfacial stresses in colloidal systems in what is arguably the most minimalist experimental configuration: the sharp interface between a colloidal suspension and its own solvent [25]. By investigating the onset of the Saffman-Taylor instability, a hydrodynamic instability ruled by the competition between viscosity and surface tension, in suspensions of PNiPAM microgels exposed to their solvent, we have shown that interfacial stresses develop at the suspension-solvent interface. These stresses will eventually decay at large enough time, because the suspension and the solvent are fully miscible, so that diffusion-driven mixing will finally erase any sharp interface. However, on short enough time scales, these stresses act as an effective interfacial tension, as it was already pointed out for the general case of miscible fluids with composition gradients more than 100 years ago, by mathematician and physicist D. Korteweg. In [25] we have shown that the effective surface tension of microgel suspensions can be rationalized in the framework of Korteweg's theory for miscible fluids. Here, we present new data on the effective surface tension between suspensions of hard, compact silica particles and their own solvent. For both kinds of particles, we also measure the bulk rheological properties, so as to compare the impact of particle structure and softness on the bulk rheology to that on interfacial stress. We find little variation of the bulk rheological properties with particle kind, while microgels and silica particles exhibit vastly different effective surface tension properties. We discuss this surprising result in view of available models for both bulk and interfacial stresses.

Sample characterization

Microgels

The soft particles are poly-N-isopropylacrylamide (PNiPAM) microgels synthesized by emulsion polymerization according to the protocol in [16] and suspended in water. PNiPAM solutions exhibit a lower critical solution temperature (LCST) close to room temperature, which results in a T dependence of the microgel size. We characterize the T -dependent particle size of our microgels by measuring the hydrodynamic radius R h with conventional dynamic light scattering (DLS) [START_REF] Berne | Dynamic Light Scattering[END_REF], using very diluted suspensions (w/w concentration c = 10 -5 ). Increasing the temperature leads to a gradual decrease of the particle size: in the temperature range 294.0 K ≤ T ≤ 309.9 K, we find that R h (T ) is well approximated by the critical-like function R h = ε(1 -T /T c) β , with T c = 304.8 K, ε = 240.85 nm and β = 0.116. All the experiments reported in this paper have been performed at T = 293.16 K, where the microgels diameter is d = 2R h = 330 nm.

We perform experiments at several particle concentrations, quantified by the effective volume fraction Φ eff = n p v, where n p is the number density of colloids and v = πd 3 /6 the volume of a single particle at infinite dilution. Note that, because microgels are soft, squeezable particles, the actual particle volume decreases at very high n p , so that samples with Φ eff > Φ rcp ≈ 0.64 or even Φ eff > 1 may be prepared, where Φ rcp is the random close packing volume fraction of hard spheres. Experimentally, only the mass fraction c of a suspension can be directly measured, by weighting a small aliquot of the sample before and after removing the solvent by evaporation. Since the microgels are highly swollen, their mass density is essentially the same as that of the solvent; consequently, c and Φ eff are proportional, Φ eff = kc. We determine the constant k using two independent methods based on the c dependence of the 3 zero-shear viscosity and of the diffusion coefficient, respectively. Viscosity measurements are performed in the range 0.0023 ≤ c ≤ 0.0980 using an Anton Paar Lovis 2000 ME microviscosimeter. Figure 1A shows that the zero shear viscosity η increases linearly with particle concentration, as predicted -for the dilute regime-by Einstein's formula:

η/η 0 = 1 + 2.5Φ eff = 1 + 2.5kc , (1)
where η 0 is the viscosity of the solvent. By fitting η/η 0 to a straight line, we determine k = 20.3 ± 0.6. Note that the viscosity data are taken for Φ eff 0.02 (see the upper axis in figure 1A), a regime dilute enough for higher order corrections to (1) to be negligible. The short-time self diffusion coefficient of the microgels, D S S , is obtained by fitting to a simple exponential decay the initial decay of the intermediate scattering function (ISF) measured by DLS, g 1 (τ ) = exp(-D S S q 2 τ ), where q is the scattering vector [START_REF] Berne | Dynamic Light Scattering[END_REF]. As pointed out in [START_REF] Pusey | Intensity fluctuation spectroscopy of charged brownian particles: the coherent scattering function[END_REF], the ISF reflects the self dynamics only if g 1 is measured at a scattering vector such that S(q, d) = 1, with S(q, d) the static structure factor. Using the value of d obtained by DLS in the very dilute limit and the Percus-Yevick approximation of S(q, d) for an equivalent hard sphere system [START_REF] Hansen | Theory of simple liquids[END_REF], we fix the scattering angle at θ = 101 deg, corresponding to q = 24.3 µm -1 , such that S(q, d) ≈ 1 at all concentrations of interest (see inset of figure 1B). Figure 1B shows that the normalized short-time self diffusion coefficient, D S S /D 0 , decreases linearly with particle concentration (here, D 0 is the zero-c extrapolation of D S S ). We model the c dependence of D S S /D 0 using the Tokuyama-Oppenheim expression [START_REF] Tokuyama | [END_REF]30]:

D S S /D 0 = 1 1 + H(Φ eff ) ≃ 1 -1.83Φ eff = 1 -1.83kc , (2)
where the approximation holds in dilute limit and H(Φ eff ) is the hydrodynamic function. By fitting D S S (c)/D 0 to a straight line, we obtain k = 20.1 ± 0.4. This value is fully compatible with that obtained by viscosimetry: the discrepancy between the two methods is less than 1%, indicating that the combination of these two low-Φ eff methods allows us to determine the absolute volume fraction of our suspensions to within 0.01-0.02 (at most), including at high Φ eff . This level of uncertainty is comparable to or better than that typically reported for colloidal suspensions [START_REF] Poon | On measuring colloidal volume fractions[END_REF]. In the following, Φ eff is calculated using k = 20.2 ± 0.5, obtained by averaging the viscosimetry and DLS values of k.

Silica Particles

Silica particles are used as a model system for hard, undeformable particles with short range, nearly hard-sphere-like interactions. The particles are Ludox-TM 50 purchased from Sigma-Aldrich and used without further purification. The particles are charge-stabilized: to partially screen long range electrostatic repulsions, we add a monovalent salt, KCl, adjusting the final salt concentration in all suspensions to 0.05 M. For very dilute suspensions, the The various symbols and colors are measurements on suspensions prepared from distinct batches of particles synthesized following the same protocol. The dashed straight lines are linear fits used to extract k, as discussed in the text. Inset of B: static structure factor S(q, d) calculated using the Percus-Yevick approximation. All DLS measurement are performed at the q vector indicated by the arrow, where S(q, d) ≈ 1 hydrodynamic diameter measured by DLS is d = 2R h = 36 nm, the same value as for saltfree suspensions. This suggests that the salt concentration is low enough to avoid particle aggregation, which may occur due to van der Waals interactions when electrostatic repulsions are only partially screened. Since rheology and surface tension measurements are performed also at high concentration, we further monitor the suspension stability in the concentrated regime (c = 0.545 w/w, corresponding to Φ eff = 0.558 as discussed below). We measure τ α , the relaxation time of the suspension at qd = 0.80, by means of multispeckle DLS, using the apparatus described in [START_REF] Masri | Aging and ultra-slow equilibration in concentrated colloidal hard spheres[END_REF], for waiting times ranging from one day up to one week after preparing the suspension. We find τ α = 4.4 ± 0.1 s independently of sample age, thus confirming that no significant aggregation occurs up to time scales larger that those of the experiments discussed here.

As for the microgels, only the mass fraction c of the silica particle suspensions can be directly measured. The two methods used to determine Φ eff and discussed above, however, cannot be applied to the silica particles, because electrostatic repulsions are not fully screened. Indeed, Einstein's law (1) holds for particles interacting solely via excluded volume interactions, while the hydrodynamic function H(Φ eff ) in ( 2) cannot be determined without a detailed knowledge of the actual interparticle potential, which is experimentally difficult to achieve. We thus use a different approach and determine Φ eff by matching the concentration dependence of the viscosity in concentrated Ludox suspensions (0.4 < Φ eff < 0.56) to that of a reference PMMA hard sphere system [33]. For solid, compact particles, one expects the effective volume fraction to be related to the mass fraction c by

Φ eff = k ′ cρ s /[ρ p -c(ρ p -ρ s )],
where ρ s and ρ p are the solvent and particle mass density, respectively. The factor k ′ accounts for deviations with respect to hard-sphere behavior, due to particle interactions: for hard spheres, Φ eff coincides with the "geometrical" volume fraction n p πd 3 /6, corresponding to k ′ = 1, while in the presence of repulsive interactions we expect the effective particle volume to be larger, yielding k ′ > 1. Figure 2 shows the result of the mapping obtained from the known values of the particle and solvent density, ρ p = 2.2 g/ml and ρ s = 1 g/ml, and using k ′ = 1.58, corresponding to a particle "interaction diameter" d 3 √ k ′ = 42 nm, about 17% larger than the hydrodynamic diameter in the limit of infinite dilution. With this choice of k ′ , a very good agreement is observed between the viscosity of hard spheres and that of the Ludox suspensions over two decades in η. To further test the robustness of this procedure, we measure the Φ eff dependence of τ α . Figure 2 shows that the concentration dependence of the microscopic relaxation time normalized by its c → 0 value is in good agreement with that of the relative viscosity, which we extend to a broader range of Φ eff using a Doolittle fit [33], η/η 0 = exp[Λφ eff /(φ eff -φ rcp )] (dotted line in Figure 2). Although the exact relation between τ α and η is still debated (see e.g. [START_REF] Poon | On measuring colloidal volume fractions[END_REF]), no strong deviations from a simple τ α ∼ η scaling are expected in the range of Φ eff probed here. Thus, the agreement between the trend of these two quantities as a function of effective volume fraction confirms the soundness of the procedure used to determine Φ eff .

Bulk rheology

Figure 3 shows the flow curves, shear stress σ vs shear rate γ, for microgels and Ludox suspensions at various Φ eff (panel A and B, respectively). The flow curves are obtained by performing steady rate rheology experiments, using a cone-plate geometry (cone diameter = 50 mm, cone angle = 0.0198 rad), except for microgel suspensions in the range 0.4 < Φ eff ≤ 1.2, for which a 25 mm-plate with a roughened surface has been used to avoid wall slip. The flow curves have been measured both by increasing sequentially the shear rate and by decreasing it, starting from its largest value. No difference were observed depending on the chosen protocol. In order to allow for a comparison with previous works, the data are presented in reduced stress and shear rate units, the former being normalized by σ T = k b T /d 3 , the typical stress created by thermal fluctuations, while the shear rate is expressed in units of the (bare) Peclet number P e = γd 2 /D 0 , where d 2 /D 0 is the characteristic time for a particle to diffuse over a distance equal to its diameter, in the Φ eff → 0 regime. At low volume fractions, for both hard and soft particles η ∼ γ throughout the whole range of shear rates investigated, indicating Newtonian behavior, as usually observed in simple fluids. As Φ eff increases, increasingly stronger deviations with respect to Newtonian flow are observed, until yield stress behavior is clearly seen, for Φ eff > 0.61 in microgel suspensions and for Φ eff > 0.59 in Ludox suspensions. Such behavior is quite general; it has predicted in recent models of complex fluids [17,34] and has been observed experimentally for a wide class of suspensions comprising both soft and hard particles [START_REF] Koumakis | [END_REF]34,35]. To better characterize the flow behavior, we fit the flow curves employing various functional forms depending on Φ eff . The fits will also be used in section 5 to obtain by extrapolation the shear-dependent viscosity in the γ regime relevant to the Hele-Shaw experiments. For Newtonian fluids, the fitting function is simply σ( γ) = η γ. For samples showing weak shear thinning, the flow curves have been fitted using a Cross-like equation: σ( γ) = η 0,s γ 1+(P γ) m , where η 0,s is the zero-shear viscosity of the suspension, 1/P is the characteristic shear rate denoting the onset of the shear thinning and m the shear thinning exponent [START_REF] Mewis | Colloidal Suspension Rheology[END_REF]. For intermediate concentrations, right below the onset of a dynamical yield stress, the shear stress is well described by a linear combination of two power laws: σ( γ) = G γg + L γl . This functional form has been successfully used to describe the flow of glassy star polymers in good solvent conditions [START_REF] Erwin | [END_REF], a model for ultra soft particles. Finally, for suspensions showing a yield stress behavior we use the Herschel-Bulkley equation [38]: σ( γ) = σ Y + λ γβ , where σ Y is the yield stress. The numerical values of the fitting parameters for all curves shown in figure 3 are reported as supplementary information [39]. Recently, there has been a lively debate on the nature of the fluid-to-solid transition in amorphous solids [START_REF] Koumakis | [END_REF]3,17,18,19,35]. In the most general scenario proposed in [3,17], up to two distinct transitions may be observed in thermal systems comprising soft particles: upon increasing Φ, a glass transition is first observed at Φ < Φ rcp , signalled by the emergence of a yield-stress plateau at low γ. The stress scale of this plateau is that of thermal stress, k B T /d 3 . Deformable particles can be further concentrated at volume fractions beyond Φ rcp , where the yield stress σ Y steeply increases to higher stress values, ruled by the energetic cost to deform a particle. Experimentally, it is typically difficult to observe these two transitions in the same system, because the range of accessible γ may be too small, or because the particle size and softness may be such that the two plateaus are not well distinct. Indeed, in most cases only one transition is seen; as a consequence, contrasting interpretations have been given to the flow curves of dense colloidal suspensions, especially in the case of soft particles such as PNiPAM microgels [17,18,19,35], where both mechanism are a priori possible. As pointed out in [3,17], it is then important to examine the numerical values of the stress and γ scales associated with the emergence of yield stress behavior. The data of figure 3 belong essentially to the "thermal sector" defined in [3,17], where P e 1 and σ Y is of the order of the thermal stress σ T . This suggest that the fluid to solid transition observed here is ruled by glassy behavior, not only for the Ludox hard particles, but also for the deformable microgels, in agreement with [35].

To compare in more detail the behavior of the two kinds of suspensions, we plot in figure 4A the Φ eff dependence of σ for three values of the Peclet number, in the regime P e << 1. Two regimes can be distinguished. For σd 3 /k b T ≤ 1 the stress shows almost no dependence on the softness of particles: the normalized stress for silica and microgels suspensions collapses on the same curve, an unambiguous evidence that both are dominated by thermal fluctuations. For σd 3 /k b T > 1, the athermal contribution to the overall stress starts emerging: The energy scale dictated by the interparticle interactions contributes to energy dissipation, producing larger stresses for the harder particles, i.e. the Ludox. Note that, because our silica particles are close to hard spheres, for the Ludox suspensions this regime extends only over a very limited range of Φ eff , close to and below Φ rcp . In the region where σ > σ T , a qualitative difference in terms of internal dynamics and bulk rheology must be expected for particles interacting via different kind of potentials. To underline this last point, figure 4B shows the yield stress extrapolated via the Herschel-Bulkley fits for the two sets of particles in the very high Φ eff regime. The yield stress of Ludox suspensions rapidly increases as a function of Φ eff on approaching Φ rcp . For the microgel suspensions, similar values of σ Y are attained, but the growth is much gentler, owing to particle deformability. We end this section on the bulk rheology by summarizing the main results: no significant differences are observed in the behavior of the two suspensions up to packing fractions close to Φ rcp , because for both systems the particles are small enough for thermal stresses to fully dominate their mechanical behavior. Only at larger Φ eff , where particle contacts are ubiquitous, does the difference in softness play a major role, as shown by the sharp growth of σ Y (Φ eff ) for the Ludox suspensions, to be contrasted to the smooth increase observed for microgels.

Viscous fingering: overview of theory

In order to measure the (transient) surface tension between the colloidal suspensions and their own solvent, we analyze the patterns formed when pushing the (less viscous) solvent in the (more viscous) suspensions. As shown in [40] for immiscible fluids, in the Hele-Shaw geometry where the fluids are confined in the thin gap between two parallel plates the interface between the displacing medium and the displaced fluid becomes instable and develops distinctive finger-like fluctuations. Crucially to our experiments, the number of fingers in this Saffman-Taylor instability is directly related to the (effective) interfacial tension. In this section, we briefly review the formalism required to retrieve the interfacial tension from the observed patterns. In the following, we will assume that no mixing occurs during the interface propagation, which is of course the case for immiscible fluids, but which also applies to miscible fluids on short enough time scales, as in our experiments. We focus on a radial geometry, where the more viscous fluid initially occupies a disk-like volume between the plates and the displacing fluid is injected through a hole in one of the plates, centered with respect to the first fluid. No lateral confinement is imposed to the fluids. At the beginning of the injection, the interface between the two fluids is stable and has a circular shape, when observed from above the two plates. As the instability develops, deviations from a circular shape are seen, which appear as "fingers" whose number and length increase with time, see figure 5a)-c). The early stages of the instability are conveniently described by decomposing in Fourier modes the fluctuations of the interface around the unperturbed, perfectly circular shape. The experimental observable is the number of fingers, which has to be related to the dominant modes in the Fourier analysis. In the last 50 years, many authors have associated the number of fingers to the Fourier mode with fastest growth rate [41,42,43,44,45,46]. This analysis is particularly convenient computationally, because an analytical expression can be derived for n f (t), the order of the mode that grows faster. However, it was pointed out in [45] that in radial Hele-Shaw experiments each Fourier mode evolves at a different time-dependent rate. Thus, the fastest-growing mode at any given time does not coincide in general with the mode that has developed the maximum amplitude. In recent theoretical and experimental studies [25,47] it has been proposed that the number of fingers observed at a given time is well approximated by the order of the mode with maximum amplitude at that time. In the following, we provide numerical evidence to support this assumption and calculate an asymptotic expression for the maximum-amplitude mode from which the interface tension may be obtained. We start from the analysis of the Saffman-Taylor instability by Miranda and Widom [45], where the perturbation around a circular interface due to the instability is decomposed in Fourier modes of (complex) amplitude ζ n (t). Assuming that the noise giving rise to the instability is a complex number ζ 0 n , with a random phase and a n-independent modulus, the time-dependent amplitude of the n-th mode can be written as [45]:

ζ n (t) = ζ 0 n K(t) (nA -1) n 2 (n -1) nA-1 exp (nA -1) 1 K(t) n(n 2 -1) nA -1 -1 . (3)
Here, A = (η 2 -η 1 )/(η 2 + η 1 ) > 0 is the viscosity contrast between the two fluids, In the very early stages, the amplitude of the instability is too small for the number of fingers to be reliably determined; accordingly, no dots are assigned to the corresponding curves ζ n (t). As soon as the fingers become visible, their number tracks very well the order of the mode with maximum amplitude, as shown by the fact that in most cases the dots correspond to the peak of the mode amplitude distributions shown in figure 5D-F, the maximum deviation being |∆n| = 1. This numerical analysis demonstrates that the number of fingers essentially coincides with n A , the order of the mode with maximum amplitude.

K(t) = [r(t)Q]/ (
To make further progress, we need an expression relating n A to the geometrical and fluid parameters, including the interfacial tension. Since an analytical expression for n A cannot be found, we derive in the following a useful expression for this mode in the limit n A >> 1. We calculate n A by solving dζ n (t)/dn = 0. Using (3), one finds that n A is the solution to

ζ 0 n K (nA -1) n 2 (n -1) nA-1 exp (nA -1) 1 K n(n 2 -1) nA -1 -1 × A( n(n 2 -1) K(nA -1) -1) + 1 K (nA -1)(3n 2 -1) -nA(n 2 -1) (nA -1) + (4) + n A n - 2(nA -1) (n 2 -1) - (nA -1) n 2 + A ln K(nA -1) n(n 2 -1) = 0
The first factor in curly brackets is non-zero for any n ≥ 1. Hence, ( 4) is satisfied only if the second factor in curly brackets vanishes, which, in the asymptotic limit n >> 1, yields

3n 2 K -3A + A ln KA n 2 = 0 . (5)
Equation ( 5) has two real solutions: n A,1 = √ KA and n A,2 = -W (-3e -3 ) KA/3, where W (x) is the Lambert function satisfying x = W (x)e W (x) . Note that n A,1 > n f while n A,2 < n f , where

n f = KA/3 (6)
is the mode with the maximum growth rate as obtained in [45]. A numerical analysis of the problem shows that the number of fingers grows with time, as confirmed by the experiments and the simulated interfaces (see e.g. figure 5B andE). Thus, at any time the mode with maximum growth rate must be larger than that with maximum amplitude. It follows that the first solution, n = n A,1 > n f , is non-physical. The final expression for the mode with the maximum amplitude is then

n A = αn f , (7)
with α = -W (-3e -3 ) ≈ 0.422. Together with ( 6) and the definitions of K and A (see (3)), (7) provides the link before the experimentally observed number of fingers and the interface tension. Before recasting this expression in a way that is more suitable to analyze the experiments, two remarks are in order. First, although equation ( 7) has been derived in the limit n >> 1, we have shown [25] that, for realistic choices of the fluid and geometrical parameters, equation ( 7) represents an excellent approximation to the full solution already for n A 2. Thus, in the next section we will identify the number of fingers observed in experiments with n A as calculated from equation (7). Second, in the above discussion we have implicitly assumed that the viscosity contrast is independent of the rate at which the invading fluid is injected, i.e. that both fluids are Newtonian. This is clearly not the case for our colloidal suspensions at large Φ eff , as seen in figure 3 (recall that η( γ) = σ/ γ). To take into account the shear-thinning behavior of the suspensions, we use the shear-rate-dependent viscosity η 2 ( γr ) for evaluating the viscosity contrast A, where γr is the shear rate at the position of the interface. This choice relies on the assumption that the wavelength of the perturbation at its onset is not drastically changed by the non-Newtonian features of the suspension. Numerical work on the Saffman-Taylor instability in a radial Hele-Shaw geometry supports this scenario [48,49], by showing that the non-Newtonian character of the fluids does not change qualitatively the instability, but just accelerates (resp., delays) its onset for shear-thinning (resp., shear thickening) fluids. This choice is also supported by previous works [50,51] on the Hele-Shaw instability between immiscible non-Newtonian fluids in a rectangular geometry, where the dynamics of the fingers was described by a generalized Darcy law where the Newtonian viscosity was replaced by the shear rate-dependent viscosity.

Having identified the viscosity of the more viscous fluid with the shear-rate-dependent viscosity of the colloidal suspensions, it is convenient to introduce a "finger function" [25] that, using ( 6), (7) and the expressions for the fluid and geometrical parameters given after (3), is shown to be proportional to the shear rate at the injection hole, γI , with a proportionality coefficient that directly yields the desired effective interfacial tension:

K * exp ≡ b r 3n 2 A α 2 -1 [4r 0 (η 2 ( γr ) -η 1 )] = 1 Γ e γI . (8)
In (8) all quantities but Γ e are known or experimentally measurable: γI = 3Q(2πr 0 b) -1 (assuming Poiseuille flow), where r 0 is the radius of the injection hole and Qb the imposed flow rate; n A is obtained by counting the fingers at the onset of the instability; r = Qbt/π is the radius of the unperturbed interface at time t; η 1 is the viscosity of the solvent, and η 2 ( γr ) = σ f it / γr , where σ f it is the fit to the flow curve discussed in section 3. Finally, the shear rate at the position r of the interface is evaluated as γr = 4r 0 γI /r, again assuming Poiseuille flow. Equation ( 8) has been shown to quantitatively capture the behavior of immiscible, Newtonian fluids [25]. It will be used in the next section to measure the effective interfacial tension between our colloidal suspensions and their own solvent.

Off-equilibrium interfacial tension: experiments

We investigate the Φ eff dependence of Γ e for both microgels and Ludox suspensions by imaging the viscous fingers in Hele-Shaw experiments. For a given Φ eff , we perform experiments at various γI , always keeping the injection rate high enough for diffusion-driven mixing between the injected solvent and the suspension to be negligible. The Hele-Shaw cell consists of two square glass plates of side L = 25 mm separated by spacers fixing the gap at b = 0.5 mm. The cell is first filled with the suspension to be studied. The less viscous fluid (the same solvent as that of the suspension, but died with 0.5% w/w of methylene blue for visualization purposes) is then injected through a hole in the center of the top plate (r 0 = 0.5 mm). The viscosity of the solvent is η 1 = 1.011 mPa s. The injected volume per unit time, Qb, is controlled via a syringe pump. Temperature is fixed at T = 293 ± 0.1 K by means of a Peltier element, with a circular hole of radius 8.5 mm for optical observation. A fast CMOS camera (Phantom v7.3 by Vision Research) run at 100 to 3000 frames s -1 is used to image the sample during injection. Typical images of the interface between the two fluids are shown in figure 6, where the instabilities are visible for both water/microgel and water/Ludox suspensions. Further images of the two-fluid interface instability at different shear rates are reported in [39]. The number n A of fingers has been evaluated by counting the number of flex points along the interface. Such number is twice the number of the fingers. To improve the accuracy of the estimate of n A , we measure it for two nearby positions of the interface: at the experimentally observable onset of the instability (corresponding to an average distance r ′ of the interface from the center of the cell) and when the average radius has grown by ∆r = 1 mm. The number of fingers is then calculated as n A = [n A (r ′ ) + n A (r ′ + ∆r)]/2; the corresponding value of the average radius used in the experimental finger function ( 8) is r = r ′ + ∆r/2.

Finally, both n A and r thus obtained are further averaged over 2 or 3 experiments performed in the same conditions (same Φ eff and injection rate). Figure 7A shows the finger function K * e vs the shear rate at the injection hole for microgel suspensions at selected volume fractions in the range 0.2 < Φ eff < 0.92. The dashed lines are linear fits to the data, showing that (8) captures very well the behavior of both Newtonian and non-Newtonian suspensions. We emphasize that strong deviations from the linear scaling K * exp ∼ γI would be observed for the high-Φ eff suspensions if the zero-shear viscosity was used to evaluate the viscosity contrast, rather than the shear-rate-dependent viscosity. Thus, figure 7A validates a posteriori the choice made in writing (8). For each Φ eff , the effective interfacial tension is obtained from a linear fit of K * exp ( γI ). Figure 7B demonstrates that very good fits are obtained for both microgel and Ludox suspensions, by showing that the law K * exp Γ e = γI predicted by ( 8) is very well verified over 4 decades. The effective interfacial tension obtained from the fits shown in figure 7 is plotted as a function of Φ eff in figure 8. The behavior of Γ e depends dramatically on the kind of particles: for microgels, Γ e increases smoothly and mildly with increasing volume fraction. For the Ludox suspensions, by contrast, the effective surface tension is smaller than for the microgels at low Φ eff (see inset), but it increases sharply as Φ eff approaches the random close packing fraction, resulting in a growth of Γ e by more than three decades over the range of packing fractions explored in our experiments. Such different behavior is also reflected in the opposite trends seen in the interface patterns, figure 6: while for microgels the number of fingers increases with Φ eff , indicating that the destabilizing effect of the growing viscosity contrast prevails on the stabilizing effect of interfacial tension, for the Ludox suspension n A decreases with packing fraction, suggesting that a massive growth of Γ e must more than offset the increase of the viscosity contrast. 8) vs the shear rate at the injection hole, γ I , for microgels suspensions with effective volume fraction in the range 0.2 < Φ eff < 0.92 (same symbols as in figure 4). The dashed lines are linear fits. B: Scaled finger function Γ e K * exp for microgel (open symbols) and Ludox (solid symbols) suspensions at various Φ eff (same symbols as in figure 4). The dashed line is the master curve Γ e K * exp = γI predicted by (8). The behavior of the microgels has been rationalized in [25] in terms of the so-called Korteweg, or square-gradient, model. We briefly recall here the main features of the model. As it was already recognized by Korteweg in 1901 [52], on time scales shorter than that of interface relaxation due to diffusive mixing, interfacial stresses between miscible fluids exist, whose effect is akin to that of an effective interfacial tension. Similarly to the theory for immiscible fluids, Korteweg's theory relates the effective surface tension Γ e to the gradient of composition across the interface. For our colloidal suspensions, one can write

Γ e = κ ∞ -∞ dΦ eff dz 2 dz ≃ κ δ Φ 2 eff , (9)
where z is the coordinate orthogonal to the interface and κ the square-gradient constant.

The last approximation holds for a linear concentration profile that increases from 0 (in the solvent) to Φ eff (the value in the bulk of the suspensions), over an interface thickness δ.

The dashed line in figure 8 shows that for the microgels (9) holds over the full range of Φ eff . This is a priori very surprising, since the square-gradient model neglects powers of dΦ eff /dz higher than two, which should become important at high Φ eff . As argued in [25], this apparent inconsistency is removed if one realizes that for the PNiPAM suspensions the relevant parameter is ϕ, the volume fraction of the monomers composing the microgels, rather than the effective volume fraction of the microgel particles. Because the microgels are highly swollen, Φ eff /ϕ ≈ 27.75, so that our experiments are actually performed in the low-ϕ limit where ( 9) is expected to be valid. As a further test of this interpretation, we have calculated in [25] the square-gradient constant κ by adapting existing polymer theories [53].

The interface thickness δ can then be obtained from κ and the prefactor in the quadratic fit to the experimental Γ e (Φ eff ) data. The value of δ thus obtained range from 364 nm to 392 nm, comparable to the average interparticle distance [25], as expected for a sharp interface, thus supporting this scenario.

For the Ludox suspension, by contrast, the quadratic law (9) fails to fit the sharp growth of Γ e vs Φ eff . This is consistent with the above picture, because the Ludox particles are compact objects, such that a sharp composition gradient is established across the solvent-suspension interface. Consequently, higher order terms not included in Korteweg's theory are likely to be important. One would still expect ( 9) to hold at low enough Φ eff ; unfortunately, the viscosity contrast between diluted Ludox suspensions and their own solvent is too low for this regime to be probed in our Hele-Shaw experiments. In the range of packing fractions of our experiments, we find that the interfacial surface tension of Ludox suspensions is well reproduced by an exponentially diverging law: Γ e = a exp(bΦ eff ), with a = 1.9 nN m and b = 35.8 (solid line in figure 8). Ongoing work in our group focuses on investigating the physical origin of this empirical law, as well as its dependence on particle size, internal structure, and interparticle interactions.

Conclusions

We have investigated both bulk and interfacial stresses in strongly driven colloidal suspensions, focusing on the analogies and differences between suspensions comprising particles with different internal structure and interparticle interactions. We find that the flow curves of soft, deformable microgels and those of suspensions of hard silica particles are very similar up to packing fractions Φ eff Φ rcp . This suggests that in our PNiPAM system the fluid-to-solid transition is driven by the same thermal glass transition as for the hard particles, and not by an athermal jamming transition. Only at packing fractions Φ eff > Φ rcp does the deformability of the microgels play a role, as revealed by the onset of a different regime where the yield stress plateau is governed by the energetic cost to deform a particle, rather than by the thermal stress.

In striking contrast to the bulk rheology, the behavior of the (transient) interfacial stress that arises when a sharp concentration gradient is imposed between the suspensions and their own solvent strongly depends on the nature of the particles. The surface tension of microgels can be rationalized in the framework of Korteweg's theory for mild interfacial concentration gradients, because for microgels the monomer concentration ϕ -and thus the interfacial concentration gradient-remains small even at high Φ eff . Moreover, a quantitative analysis [25] inspired by work on polymer solutions [53] allows one to determine the square-gradient constant κ, and therefore to account for the absolute magnitude of Γ e . This result is based on the calculation of the free energy cost associated to the reduction of the degrees of freedom of the polymers chains near the interface, due to the additional constraint imposed by the presence of a concentration gradient [53]. Thus, both the applicability of Korteweg's theory to PNiPAM suspensions and the absolute magnitude of Γ e are ultimately the consequence of the polymeric nature of the microgels, rather than their deformability per se. It would be interesting to further test these ideas by measuring Γ e for suspensions of soft but nonpolymeric particles, e.g. using emulsions. Our understanding of interfacial stresses for hard particles is much less advanced. We find that Γ e increases exponentially with Φ eff , as opposed to the quadratic growth seen for the microgels and predicted by Korteweg's theory. Deviations from the square-gradient law are not unexpected for systems exhibiting strong concentration gradients, as is the case for the sharp interface between the solvent and a dense suspension of compact particles. However, the origin of the precise functional form of Γ e , its magnitude, the influence of the particle material, softness, and interparticle interactions are still open question that we are currently investigating.

Figure 1 .

 1 Figure 1. Relative viscosity (panel A) and normalized short-time self diffusion coefficient (panel B) of PNiPAM microgel suspensions at different concentrations c (w/w).The various symbols and colors are measurements on suspensions prepared from distinct batches of particles synthesized following the same protocol. The dashed straight lines are linear fits used to extract k, as discussed in the text. Inset of B: static structure factor S(q, d) calculated using the Percus-Yevick approximation. All DLS measurement are performed at the q vector indicated by the arrow, where S(q, d) ≈ 1

Figure 2 .

 2 Figure 2. Relative viscosity (red full squares) and normalized microscopic relaxation time measured by DLS (red empty squares) of Ludox suspensions as a function of the effective volume fraction Φ eff . Reference hard-sphere data for the relative viscosity have been taken from [33] (black full squares) and fitted by the Doolittle equation (dotted line).

Figure 3 .

 3 Figure 3. Flow curves: stress versus shear rate under steady shear for samples with different effective volume fractions. MICROGELS (Panel A): from bottom to top, Φ eff = 0.20, 0.27, 0.33, 0.46, 0.55, 0.61, 0.92, and 1.2. LUDOX PARTICLES (Panel B): from bottom to top, Φ eff = 0.41, 0.46, 0.540, 0.558, 0.563, 0.573, 0.591, 0.607, 0.621, and 0.632.

Figure 4 .

 4 Figure 4. A: Normalized shear stress as a function of Φ eff for three different Peclet numbers in the thermal sector (P e < 1), for microgels (empty symbols) and solid particles (solid symbols). Some of the values reported here have been obtained by extrapolation using the fitting functions described in the main text. B: Φ eff dependence of the normalized yield stress for microgels and Ludox suspensions obtained via Herschel-Bulkley fits.

  2πβ), where r(t) is the distance from the center of the cell of the unperturbed fluid-fluid interface, Q is the area invaded by the injected fluid per unit time, and β = b 2 Γ/[12(η 1 + η 2 )], with b the cell gap and Γ the interfacial tension between the two fluids. Note that (3) only holds for nA > 1. In our experiments, the viscosity contrast is sufficiently large for this inequality to be fulfilled for n ≥ 2. We show in the top row of figure 5 a top view of the theoretical interfaces calculated using equation (3) for three sets of realistic fluid and geometrical parameters. In each panel, the interface is plotted at several times after starting the injection of the less viscous fluid, from t = 10 -4 s in steps of 0.1 s. The bottom row shows the corresponding temporal evolution of the amplitude of the Fourier modes, as a function of the mode number. The dots in the bottom-row plots indicate the number of fingers observed in the patterns shown in the top row.

Figure 5 .

 5 Figure 5. Theoretical time evolutions of the fingering patterns (top row) and correspondent Fourier mode amplitudes (bottom row). The time between successive data sets is 0.1 s. The dots in the bottom row plots indicate the number of fingers visible in the patterns in the top row. The interface patterns are generated using (3), with the following sets of parameters: A,D: Γ = 4 mN/m, V = 8 • 10 -9 m 3 /s, η 1 = 10 -3 P a • s, η 2 = 1 P a • s, |ζ n (0)| = 10 -8 m. B,E: Γ = 4 mN/m, V = 8 • 10 -9 m 3 /s, η 1 = 10 -3 P a • s, η 2 = 4 P a • s, |ζ n (0)| = 10 -9 m. C,F: Γ = 4 mN/m, V = 8 • 10 -9 m 3 /s, η 1 = 10 -3 P a • s, η 2 = 6 P a • s, |ζ n (0)| = 10 -9 m.

Figure 6 .

 6 Figure 6. Solvent-suspension interface for PNiPAM microgels and Ludox suspensions, for various effective volume fractions and shear rate at the injection point as indicated by the labels.

1 ]Figure 7 .

 17 Figure 7. A: Finger function K *exp defined in(8) vs the shear rate at the injection hole, γ I , for microgels suspensions with effective volume fraction in the range 0.2 < Φ eff < 0.92 (same symbols as in figure4). The dashed lines are linear fits. B: Scaled finger function Γ e K * exp for microgel (open symbols) and Ludox (solid symbols) suspensions at various Φ eff (same symbols as in figure4). The dashed line is the master curve Γ e K * exp = γI predicted by(8).

Figure 8 .

 8 Figure 8. Effective interfacial tension Γ e between colloidal suspensions and their own solvent as a function of colloid effective volume fraction. Open symbols: microgels (data from[25]); solid symbols: Ludox suspensions. The dashed (solid) lines are quadratic (exponential) fits, as discussed in the text. Inset: same data plotted in semilogarithmic scale.

Acknowledgments. We acknowledge fruitful discussions with J.-L. Barrat and L. Berthier. This project was supported by ANR under Contract No. ANR-2010-BLAN-0402-1.


P Cuenot1

E Jenn1

E Faure2

N Broueilh2

E Rouland1

An Experiment on Exploiting Virtual Platforms for the Development of Embedded Equipments

Keywords: Virtual platform, SystemC, Verification. Simulation, Heterogeneous environment

Virtual engineering methods and tools based on simulation have become a privileged mean to reduce time-to-market and product cost. However, design and verification activities still need to be improved to manage the ever increasing complexity of electronic products and their interactions with heterogeneous environments. In particular, an important challenge is to master the real time properties of the product composed of interacting hardware and software components.

In this paper we propose a pragmatic approach to use virtual platforms to verify gradually and accurately the properties of a system under design. We illustrate the approach on an example.

Introduction and Motivation

New embedded electronic systems impose a new leap in product integration. The progress of the System on Chips (SoC) market share clearly illustrates this need. In parallel, embedded electronic systems or Electronic Control Units (ECU) need to be optimized to integrate new functions usually developed using analog parts at an acceptable cost. Basically, the challenge is to maximize the usage of SoCs while keeping reasonnable development costs. Embedded electronic systems are usually reactive: they are part of some control loop closed through their environment. The temporal properties of the elements of this loop are of prime importance as they determine the overall performance of the implemented function (response time, stability, etc.). Accordingly, those properties must be carefully monitored and the product (including drivers, IPs, sensors, actuators…) shall be carefully designed to comply with those properties. In this context, being able to simulate a complete control chain within its environment is of major interest, and the benefits of the approach are as high as it can be applied soon, continuously and "smoothly" during the development process. Towards this goal, there is a strong need to move from segregated hardware and software simulation to system-level simulation. This requires the modelling and simulation of the system's components and environment at various levels of abstraction, representativeness, and accuracy. The SystemC language and simulation kernel [START_REF]SystemC core langage definition, Accelera Systems initiative standards[END_REF] provide such capabilities. SystemC is supported by an open source ecosystem organized in the context of Accelera Systems Initiative [START_REF]SystemC core langage definition, Accelera Systems initiative standards[END_REF]. Recently, the Electronic Design Automation tools industry (EDA) extended commercial offers by integrating SystemC. The tool environment offers modelling services on the top the language, scripting capabilities to automate test execution, co-simulation interface to interact with external tools or actual electronic prototypes, etc. Despite (or because of) the large offer of technical methods and tools, some methodological guidance is required to ensure a safe, high-quality, and cost-effective development and verification process. Therefore, we (i) propose an iterative process to optimize the use of SystemC, (ii) apply it on part of a system, and (iii) show its benefits. Focus is placed on verification activities carried out with the system's environment to demonstrate the preservation of the components properties during the successive development phases. Our paper is structured as follows. Section 2 introduces the virtual platform concept and gives an overview of some significant related works. Section 3 presents the proposed approach for an iterative development and verification process. In section 4, the experimental setup is explained and section 5 gives some evaluation results. Finally, the conclusion reports feedback on methods and tools evaluation by application of the process and draws perspective for future work.

Related Work

A virtual platform is a hardware simulator executing embedded software. The hardware simulator is usually built on top of an Instruction Set Simulator (ISS) of the processor connected via buses to memory and peripherals such as timers, general I/O, communication interface, etc. The ISS may be implemented in SystemC or in any other general purpose language. In a typical configuration, communications are modelled using SystemC Transactional Level Modelling (TLM) in order to achieve high simulation speed. Peripherals and memory are register accurate, they communicate according to the TLM paradigm, and their behaviour is implemented in SystemC/C++. The application software is the actual code compiled for the target processor. The SystemC non-preemptive simulation kernel orchestrates the execution of all components of the platform. Among typical examples of virtual platforms, we can mention the Infineon Tricore TM SoC [START_REF] Koppelmann | An open and Fast Virtual Platform for Tricore TM based SoC Using QEMU[END_REF] based on the QEMU ISS or the SockROCKET LEON3 virtual platform developed by ESA [START_REF]ESA SockROCKET virtual platform[END_REF] implemented in Python with SystemC/TLM buses and peripherals. Note that the ISS can be developed by silicon suppliers from proprietary architecture modelling languages, such as Freescale's ADL/uADL [START_REF] Freescale | uADL open source project[END_REF], and then be integrated with peripherals to build a complete platform. Virtual platforms have been experimented on various industrial use cases, such as automotive power train applications [START_REF] Cuenot | Embedded Software V&V using virtual platforms for Powertrain application[END_REF]. Those experiments have demonstrated the need for appropriate methods during the development of the platform components and their integration.

In the industry, the different uses of virtual platforms map to the organisation of the supply chain: design and verification of SoC on the silicon suppliers side (component), design and verification of software on the equipment suppliers side (system). On the system side, virtual platform are used to estimate the performances of hardware and software architectures, and perform early verification activities. In particular, debug and verification phases are simplified thanks to the high observability and controllability of virtual platforms compared to real hardware. Furthermore, using virtual platforms moves the simulated hardware parts out of the critical path: hardware development phases can start once the definition of the hardware has matured and has been (virtually) validated. Finally, virtual platforms also provide: -A high configurability and flexibility allowing a particular platform to be configured and elaborated within minutes provided that the models are available. -A capability to integrate models with heterogeneous abstraction levels (in particular temporal abstractions thanks to the versatility of Sys-temC/SystemC-TLM/ SystemC-AMS) In its Return On Investment (ROI) analysis on electronic system level design, Synopsys [START_REF] Schirrmeister | System Level Market trends[END_REF], one of the main EDA tool supplier, announces cost reduction opportunities by reducting the number silicon iterations and a productivity increase of x2-x5 for software development.

To reach such a ROI, the virtual platform design flow shall be optimized so as to (i) minimize the cost of models and (ii) maximize the credit one may take from verifications performed using those models. An approach consists to use simulation models as a form of an "executable design artifact" that is refined all along the development process, from the implementation agnostic logical level down to the physical software and hardware levels. Obviously, this approach is conditionned by the quality of the model and particularly on the properties preserved by the model. Those qualities must be clearly stated and become part of a contract binding the provider and the user of the model. A specific care shall be taken on timing properties of the hardware model and of the environment impacting the overall system behaviour. The Socket project [START_REF] Lefftz | A design flow for critical embedded systems[END_REF] has proposed such a design flow for the development of critical embedded SoCs. The development steps maps to the SystemC/TLM modelling styles. This allows to co-simulate hardware and software, and to introduce and verify properties gradually (in particular temporal properties). The top level design of the SoC architecture is focused on hardware bus traffic optimization, not on the complete set of properties allocated to the SoC. In this project, the modelling of the environment is limited since the systems considered are not reactive. The COMPLEX project [START_REF] Herrera | The COMPLEX methodology for UML/MARTE Modeling and design space exploration of embedded systems[END_REF] uses UML/MARTE to model and explore the design space of embedded systems. Power and performance are the exploration criteria. The various abstractions of the processor micro-architecture, of the SoC internal buses, etc. allow an early and efficient simulation. Those abstractions also impair accuracy and raise sensitivity problems (in particular with target compiler optimization versus native compiler and internal SoC bus communication control). Compared to Socket, COMPLEX placed focus at system level where interconnected ECUs communicate through communication buses. The precision in low-level modelling is not considered besides the use of existing SystemC component libraries. Additionally the verification activities are not formalized. With respect to the previous works, our approach is aimed at covering a larger modelling spectrum, from high level models down to behavioural hardware models. Emphasis is placed on verification of realtime properties considering the ECU's environment.

An Approach for the Development and Use of Virtual Platforms

To benefit the virtual platform's "good properties", the objective is to (i) minimize the development cost of the virtual platform and (ii) maximize the usage of the virtual platform.

To achieve the first objective, one shall (i.a) maximize the reuse of existing models (possibly reusing them from previous developments), (i.b) maximize the automatic generation of platform models (or skeletons) from existing design models, (i.c) develop simple models covering the necessary and sufficient features and details with respect to the validation and verification objectives. These topics will be addressed in Section 4.d.

To achieve the second objective, the strategy is twopronged: (ii.a) the "most expensive" problems shall be identified and tackled first, (ii.b) the model shall be detailed up to the point where the necessary and sufficient precision and accuracy are obtained to solve (ii.a).

Here, we propose an informal approach somewhat similar to a Failure Mode and Effect Criticality Analysis (FMECA) applied on a development process. This approach takes into account: the cost of evaluations including the cost of model development, the cost of model execution and the risk inherent to a "sloppy" evaluation. According to this interpretation, the frequency of occurrence of an error is related to the quality of the measures (its accuracy and precision). The gravity integrates the potential impact of the error on the development process (e.g., the number and nature of the activities to be redone), and the impact on the product under design. Once the analysis is complete, reducing the overall risk consists to reduce the probability of occurrence of the evaluation error (e.g., by refining the model in order to account for phenomena that have a significant impact on the behaviour of the system), (ii) reduce the gravity of the error by improving the robustness of the process (e.g., by introducing margins), (iii) detect evaluation error by adding additional checks. In this paper, focus is placed on the temporal properties, so the trade-off between cost and accuracy/precision is essentially focused on the modelling of temporal aspects. We consider the three levelsor programming styles -defined by SystemC-TLM: untimed (U or functional), loosely-timed (LT), and approximately-timed (AT). As the design improves, the verification objectives get more and more focused and may eventually require a fine grain temporal modelling. This change is determined by the property to be verified.

Failure modes

Contrary to hardware devices, there is no common, standardized list of failures at functional level. Candidate failure modes are identified and selected on the basis of expertize, experience, etc. Here, focus is placed on temporal errors such as under and over estimation of delays, non-compliance with SW/HW interaction sequencing or hardware design constrainst. The objective is to determine which verification technique to apply by considering the design faults, their effect on the system and the cost of the detection / mitigation means.

Criticality

Our estimation of criticality is (roughly) estimated by a cost. The principle is trivial: the cost of the verification means shall be somehow commensurate or related to the direct and indirect costs of the error. The cost of the error is roughly estimated by its criticality, which depends on its probability of occurrence and the costs of its effects including the cost of detection means (technical and human) and the cost of error elimination (redesign and refactoring,…).

Error propagation

With respect to a classical FMECA, the approach differs because design errors propagate across design artefacts (space) and across design phases (time). In particular, an error in the dimensioning of a parameter in phase 𝑖 may impact the design choices done in phase 𝑖 + 1, and so on. Any verification approach is basically aimed at preventing such propagation by (i) detecting design errors as soon as possible in the design process and (ii) as soon as possible in the dependency chain that relate design artefacts. This "approach" has been applied on some simple functions of our experimental setup. The corresponding use cases are described in section 4.b and 4.c. Back to the important question of the ROI of such multi-viewpoint, multi-level approach, the following questions shall be answered: -What is the cost of developing (all) those additional models? The answer obviously depends on their complexity. Hence, the development effort for a SystemC-TLM timed model of a simple timer is lower than one men×month, while a complex ISS can require more than twenty men×months. -If models are develop by some third party, how do we formalize the "contract" that binds the model provider to the model user? Stated differently, how can we specify the domain to obtain significant results with the virtual platform model? -Finally, how does this additional cost compares to the gain due to the early validation? In the sequel of the paper, we propose to answer those questions by applying the virtual platform approach on a small example: an autonomous rover.

The experimental setup

In order to evaluate the proposed development process and supporting tools, and estimate its actual benefits, we use it for the development of a small mobile vehicule, or "rover", named "TwIRTee".

a. The TwIRTee demonstrator

TwIRTee is a three-wheeled autonomous rover fitted with a camera and various other sensors (odometry, positioning,…). Its operational role is very simple: (i) move itself on some predefined tracks from a point A to a point B (a "mission") while avoiding other rovers.

TwIRTee is developed within the INGEQUIP project at the Toulouse Institut de Recherche Technologique (IRT) Saint-Exupéry. IRTs are new research structures established under the auspices of the French Agence Nationale de la Recherche (ANR) aimed at favouring the transfer of innovation from laboratories to industries. TwIRTee is designed so as (i) to cover the major topics addressed in the project namely: early validation, architecture exploration, performance prediction, and formal verification. Furthermore, it is aimed at covering issues, functional and architectural elements specific to three industrial domains: automotive, space, aeronautics. Accordingly, the missions, functions and the architectural elements are determined so as to tackle or exercise one or several issues: for instance, the "localization" function relies partially on imaging so as to exercise hardware / software space exploration and codesign; the highly redundant architecture provides the experimental setup to perform early performance evaluations (including dependability). An overview of the TwIRTee plaftorm is given on Figure 2: the computing platform is composed of 2 COM/MON channels that host the main "mission" functions and one channel dedicated to power supply generation and motor control. A clock synchronization protocol is implemented and distributed on each ECU communicating by the CAN network. The rover displacements is achieved by the motor control ECU controlling a two-wheeled powertrain composed of 2 CC motors, 2 reduction gearboxes, 2 quadrature encoders, 2 wheels. The setpoint for motor regulation is selected from the 2 commands (COM) and monitoring (MON) channels. The methods introduced in Section 3 are applied on two simple functions: the clock synchronization (Fck) and the PWM motor control (FPWM).

b. Clock synchronization

The clock synchronization protocol is used to provide the rover's computation units with a "common" time reference. The protocol has been proposed in [START_REF] Rodrigues | Faulttolerant clock synchronization in CAN[END_REF]. It is built on top of the CAN network.

Principles

The common time referenceor virtual clock (𝑣 ) satisfies the precision, rate and accuracy properties. The precision property states that no two synchronized virtual clocks may differ by more than a given value. For instance, at any time 𝑡, any virtual clock shall show a time "less than 100𝜇s" from any other virtual clock. More formally, if nodes 𝑘 and 𝑙 participate to the protocole:

∃𝛿 𝑣 : |𝑣 𝑘 (𝑡) -𝑣 𝑙 (𝑡)| ≤ 𝛿 𝑣 (P1)
In practice, the achievable precision depends on two main parameters Γ tight and Γ max agree : where -𝑇 is the resynchronization period, 𝜌 is the physical clock drift -Γ tight is the network propagation delay (around 10𝜇s), Γ max agree is the agreement delay which depends in particular on the number of tolerated faults and on the background traffic of higher priority. First, let's analyse the failure modes, their "probabilities" of occurrence, and their effects. Failures -FF1M1: Erroneous Γ max 𝑎𝑔𝑟𝑒𝑒 , underestimation -FF1M2: Erroneous Γ max 𝑎𝑔𝑟𝑒𝑒 , overestimation -FF2M1: Erroneous Γ tight , underestimation -FF2M2: Erroneous Γ tight , overestimation Probabilities of occurrence -OF1M1: VERY HIGH, because (i) many factors contribute to the communication times, and (ii) the dynamic behaviour of CAN is not trivial (see [START_REF] Broster | Timely use of the CAN protocol in critical hard real-time systems with faults[END_REF] and [START_REF]Road Vehicles -Controller area network (CAN) -Part 1: data link layer and physical signalling[END_REF]). -OF1M2: LOW, because the mode for fault F1 is much likely F1M1. -OF2M1: VERY LOW because this parameter is part of the specification. -OF2M2: VERY LOW (same reason as F2M1). For space reason, we do not consider failure modes F2 any longer. Effects and cost impact -EF1M1: Actual accuracy lower than expected. Cost impact is MODERATE: (i) the system being redundant (MASTER/SLAVE, COM/MON) many functions rely on the synchronization precision in particular via discrepancy margins, confirmation times are affected. However, the risk for a large effect on the synchronization is low thanks to the 𝜌 1 factor in (P2). -EF1M2: Actual accuracy greater than expected.

𝛿 𝑣 = 𝛿 𝑣 + 2𝜌𝑇 (
Cost impact is LOW: the network and CPU loads are higher than necessary but no rework is expected. (Note that, generally speaking, the impact could be very high because it could lead to select and oversized computing platform and / or network. In our very case, there is no risk of such effect). In the absence of dedicated detection means, FF1M1 and FF1M2 are likely to be detected only in operations. Probability of detection -DF1M1: VERY LOW because Γ max is a worst-case situation that is difficult to observe in operation. -DF1M2: VERY LOW because the effects are hardly visible.

Cost of correction -CF1M1: VERY HIGH -CF1M2: VERY HIGH

1 𝜌 is in the range of 10 -6 s/s From the combination of a MODERATE cost, a VERY LOW detection probability and a VERY HIGH detection / correction cost, we decided that it was worth adding a new detection means with a MODERATE cost and HIGH detection coverage. Consequently, we decided to test the clock synchronization protocol on a bit-level virtual model of the CAN capable of simulating the actual effects of background traffic and error occurrences. (Note: we consider that the behaviour of the network in the presence of fault is simulated. However, in the particular case of CAN, analytical models of the CAN bus latency are already available (e.g., [START_REF] Broster | Timely use of the CAN protocol in critical hard real-time systems with faults[END_REF]), but simulation models allow to consider faults models of arbitrary complexity.

c. PWM motor control

Power is delivered to the rover motors via a H-bridge (4 transistors, see Figure 3) controlled by a PWM signal. The PWM duty cycle is computed by the motor regulator from the speed set point provided by the rover guidance controller and the actual speed measured using the wheels' optical encoders. The wheels' optical encoders generates quadrature signals acquired as frequencial input and then integrated for speed determination. The PWM controller was assessed under two perspectives: PWM timing and PWM design constraints both having effect on real time properties.

PWM control generation

The PWM generation device hold the following functional requirements (P4 -non exhaustive list) -Frequency and resolution step of PWM control: 10kHz frequency with a 0.1 µs step) -Latency for application of new PWM value (next cycle) -Immediate desactivation of active state of the PWM (lower than 1 ms) Compliance to the previous requirements may be achieved in many different ways, including: -Pure software implementation toggling an output port (resolution will be certainly an issue) -Pure hardware specific implementation with single register interface and fixed frequency (too specific solution) -Mix mode with a simple timer, control by software interrupt and latch of an output (resource will be an issue) -Hardware dominant solution with multiple register interface but "reasonable" or limited resource (most realistic in this basic example) Attached to selected scenario we have also a list of design constrainst (C4) stemming from domain experience, available resources, etc. For instance : -Register shall be 32-bit wide. At most 4 registers must be necessary (duty cycle value, frequency value, cycle counter and register control) -The timer frequency range shall be in [100Hz,100KHz] (10kHz for motor control) -The effect of the frequency and duty cycle change shall only occur at the next PWM cycle For the evaluation of the motor controller we will consider the failure mode "wrong design" (F2M1). The interaction between the PWM and the H-bridge (the device that physically controls the power signal) raises another possible failure mode. As the H-bridge provides no overcurrent protection, care shall be taken not to switch transistors located on the same side of the bridge on "at the same time" in order to avoid shortcuts. This situation can be prevented by introducing a silence time between the commutations of opposite transistors. The duration shall in particular take into account the transistor switching time. The respective failure mode are noted F1Mx. The following property shall hold:

∀ 𝑡𝑏 , 𝑡𝑓 𝑗 ∶ |𝑡𝑏 -𝑡𝑓 𝑗 | ≥ Δ𝑡 𝑑𝑧 (P5)
Where 𝑡𝑏 (resp. 𝑡𝑓 ) be a time at which signal "backward" (resp "forward") is active, and Δ𝑡 𝑑𝑧 is the duration of the "silent zone" where no transistor is active. Failures modes -FF1M1: Erroneous Δ𝑡 𝑑𝑧 , underestimation -FF1M2: Erroneous Δ𝑡 𝑑𝑧 , overerestimation -FF2M1 : Wrong design Probability of occurrence -OF1M1: LOW -OF1M2: LOW -OF2M1: MEDIUM (component are selected to fulfill all forecoming applications) Effects -EF1M1: During a rotation sense change, the opposite MOSFETs of the H-bridge transistors are on at "the same time". Cost impact is VERY HIGH: this configuration is basically a shortcut. Depending on the duration of the shortcut, the dissipated energy may leads to the destruction of the two transistors and the power supply. No other function is affected by the error.

-EF1M2: During a rotation sense change, the PWM signal is delivered slightly later to the H-bridge.

Cost impact is NEGLIGIBLE: the rotation of the wheel is slightly delayed which has no further impact on the rover's capabilities. -EF2M1: In case of wrong design with non capable component the function performance must be downgraded. In the worst case, a complete redesign may become necessary. Cost impact is HIGH because the problem will be detected during verification but it will have an impact on overall product planning. In the absence of dedicated detection means, FF1M1 and FF1M2 are likely to be detected only in operations. FF2M1 will be detected lately during product verification meaning high effort for redesign. Probability of detection -DF1M1: LOW. Depending on the duration of the "short-circuit", the number of forward/backward commutations, the error may stay undetected for a long time. However, it may eventually reduce drastically the lifetime of the transistor / power supply. -DF1M2: VERY LOW. In practice, we created a virtual platform to host the PWM driver software and implemented a system-C observer to check (P5).The virtual platform shall provide a representation of time "compatible" with the property at stake. Here, we used a SystemC AT model of the PWM hardware driver. For the "wrong design" fault model, the combination of HIGH cost, HIGH detection and HIGH (to VERY HIGH) correction cost leads to the introduction of a new detection means aimed at securing later development phases. As the functional requirement (P4) and design constrainst (C4) are expressed in terms of hardware and software interactions, the virtual platform shall at least provide a LT abstraction. To facilitate the verification of the functional properties (P4) a functional model is developed (initiating also the test bench environment). From this level, requirements and associated properties are refined and al-located. During the decomposition process, some requirements may become "contracts" binding the different components of the motor control chain. Figure 4 depicts some of these contracts: -The PWM controller shall ensure a silence-period greater than Δ𝑡 𝑑𝑧 and the H-bridge transistors shall switch is less than (e.g.,) Δ𝑡 𝑑𝑧 /100 (P5). -The power dissipated by the electrical motor shall be less than 15W. Therefore, the PWM duty cycle shall never be greater a given ratio for a given time. This contract propagated from motor power dissipation constrainst is not considered in the article.

d. The virtual platform

VLAB™ 2 is used to develop the models and build virtual platforms. It provides a programmable and interactive environment for the assembly, configuration, programming and operation of electronic system level simulations, such as virtual system prototypes and virtual platforms, as well as other applications (see Figure 5) . A virtual platform integrates together simulation models and other simulation objects, scripts, tools, test and infrastructure software, and target software. The tool framework facilitates the creation of models thanks to a Genesis library compliant with SystemC and IP-XACT standards. It allows to capture model structure and connectivity and to automatically create C++ skeleton for the implementation of the behavioral part of the model. This library is also available in Python meaning models can be developed and debugged in a much faster and simpler way than in C++.

Figure 6: Genesis Framework

Finally, the tool provides an integrated and interactive execution platform leveraging again Python capabilities. It provides a simple yet efficient user experience for assembling virtual platform, debugging virtual platform (setting HW & SW breakpoints) and scripting test scenario (including fault injections in a non intruisive way) .

Modelling and assembling a virtual platform in Python are the key functionality used to support the implementation of the case study.

Implementation approach a. Clock synchronization

We used ASTC VLABTM "CAN toolbox" to create a model of our computation and communication infrastructure (5 nodes) and to inject faults during the execution of the clock synchronization protocol. A CAN node comes with two models, at token and bit levels.

The toolbox provides an API for the configuration and control of nodes, including activation, frame transmission and frame reception. Python scripting allows to access directly to API, to send CAN frame and control the bit engine in order to introduce error in the CAN frame. These features have been used to exercice the clock synchronization protocol in situations where bitlevel errors lead to multiple retransmissions of the same message. Such situations would have been difficult to obtain using the actual hardware.

b. The Motor control

The PWM controller scenario, as elaborated in section 4.c, is realized using VLAB TM . Some parts of the complete model were developped using SystemC and Python modeling capabilities while some other parts were directly taken from the hardware library (for TLM transactor or MPC5674F AT models). The modelling scope for contract verification of the PWM controller is enlarged to the overall motor controller feature. For the sake of the demonstration, we consider that the C code of the motor controller was generated from the same Simulink® model used to design and validate the controller.

The development phase of the motor controller corresponds to the three predefined modelling style U, LT and AT.

Untimed model

The Untimed model is a functional model. It describes how the PWM signal is generated from the rotation speed setpoint and the actual speed measured from the optical encoders (see Figure 7 below). The Hardware Abstraction Layer API (HAL) is defined at this level. It allows to implement a static interface between application layer and driver abstraction using physical data interface with the motor regulation algorithm (% for PWM and speed in km/h for integration of quadrature signal).

Motor Control

Python Comp.

Driver Abstraction

Freq. In This first level of modelling is a necessary step to build the subsequent models. We take benefit of it to perform a first verification of the design with respect to properties P4 and P5.

PWM Motor Regulation

To do so, a test environment (or testbench) is built as shown on Figure 8. In particular, property P5 (dead zone) is checked using a dedicated SystemC/python component (P5 checker on the Figure 8) that is connected to PWM output signals via a monitor component). The test driver generates a scenario where wheels are moved forward and backward.

To close the control loop, a very simple model of the motor / encoder devices is developed and integrated in the platform. In this model, the frequency of the quadrature signals is considered directly proportional to the value of the PWM duty cycle (this is a good approximation as far as the frequency of the PWM signal is sufficiently high). This model is used to get a first insight on the performances of the regulation (response time, stability). In the future we will investigate co-simulation with Simulink® dynamic of motor itself. The HSI is implemented using five 32 bits registers:

VLAB TM

System Under Test

-CNT running counter register (range of 24 bits) with 10ns resolution (range from 6Hz to 100 MHz, capable for 10KHz with 0.1% of precision) -A and B compare registers for PWM control (range of 24 bits) -CCR control register with enable bit (UCPREN bit 6), polarity bit (EDPOL bit 24) and fixed buffered mode (OPWFMB matching eMIOS definition bit 25-31) -CSR status register with overrun bit (FLAG bit 31) -Sequence transfer for buffered register A,B at the end of the PWM cycle and guarantee of B always greater than A. The objective is not to implement the complexity of an hardware IP like the eMIOS, but build a simplified IP fulfilling the fonctionnality with respect of a reduced HSI accessed via TLM2 access (untimed). The motor control is now split in two python components as depicted in Figure 9. This split allows the observability using VCD. The python components are then instanciated with the test bench architecture as in Figure 8. For the testbench, the interface for the motor control is identical as the functional model. The P5 checker can also be reused to demonstrate the dead zone contract. See Figure 11 for trace of the results.

The HSI interface contract can be demonstrated thanks to respect of the API resolved during virtual platform building and sequence demonstrated by instrumentation capability to trace TLM transaction (register access and value transported).

Approximated Timed model

The last model integrates approximated timing on bus communication and on hardware IP resource component access (internal definition and average processing time in hardware IP). Model can also complete the register interface when extra control registers are necessary for pre-implementation constraint, such as merging several LT models into a more complex and configurable hardware component. In any case the predefined HSI contracts shall be preserved: at least the same registers are used to interact with the component, according to the same sequence (protocol). A use case more relevant than timer IP could be an image processing hardware accelerator merging LT models with modeling of internal timing interconnection.

In our experiment, the LT timer IP is mapped to the MPC5674F eMIOS AT model. The Timer IP is connected to a bridge and to the memory bus of the Pow-erPC core via the peripheral bus (see Figure 10). The target software is the binary software, including final drivers, application software compiled for the µC target.

Motor Control

Conclusions and future work

In current SystemC ecosystem, the adequacy of model representativity for verification objective and the quality of the simulation models curbs the large deployement of virtual platform. Our approach provides a first level of guidance and formalization of the modeling process.

The "FMECA-like" approach introducesto some limited extend, howeverobjective criteria to determine whether a refined model is necessary or not. It allows to start the modelling phase at the adequate abstraction level for the verification of the selected property.

The formalization of the development and the test of an architecture compliant with SystemC modelling standard allows a sound design covering hardware software codesign process. We demonstrated that test elements can be reused all over the model refinement. Moreover, the contract approach provides backbone for decomposition and guarantee on the coverage of the requirements. The tool automation features simplify the platform models generation for component interfaces description and component bindings. However, for the creation of a new component model, the remaining manual operations (code implementation, platform integration and test campaign) still represent around 90% of the development costs. In actual practice the elaboration of simulation models follows a typical topdown approach based on successive abstractions. Such refinement has been formalized through the control motor modelling experiment. More important, the models have been continuously verified. In our case, the development and testing of two additional abstracted models is only estimated at 10 to 20% of the overall effort. Balanced to the model development effort, the proposed methodology with virtual platform provides strong benefits on the overall product development.

The contract and the model based communication facilitate inter domain communication and improve the design quality of the product. It helps to reduce significantly the re-work iterations. The effort gained can be estimated at 10 to 20%. This was demonstrated in our experimentation and also on porting the OS Trampoline [START_REF]Trampoline OpenSource RTOS projet[END_REF] on MPC5674F microncontroller for twIRtee. The early hardware software integration by simulation reduces operational set up and validation. The approach proposed in this paper will be experimented for the development of other parts of the INGEQUIP project demonstrator. In particular, the camera-based line tracking function will be developed according to the same approach. This function requires a lot of computing power, focus will be placed on the issue of hardware/software design space exploration.

  detection : cost of detection / correction

Figure 1 :

 1 Figure 1: FMECA overviewThis "approach" has been applied on some simple functions of our experimental setup. The corresponding use cases are described in section 4.b and 4.c. Back to the important question of the ROI of such multi-viewpoint, multi-level approach, the following questions shall be answered: -What is the cost of developing (all) those additional models? The answer obviously depends on their complexity. Hence, the development effort for a SystemC-TLM timed model of a simple timer is lower than one men×month, while a complex ISS can require more than twenty men×months. -If models are develop by some third party, how do we formalize the "contract" that binds the model provider to the model user? Stated differently, how can we specify the domain to obtain significant results with the virtual platform model? -Finally, how does this additional cost compares to the gain due to the early validation? In the sequel of the paper, we propose to answer those questions by applying the virtual platform approach on a small example: an autonomous rover.

Figure 2 :

 2 Figure 2: Overview of the TwIRTee equipment

Figure 3 :

 3 Figure 3: Motor H-Bridge control

  -DF2M1: HIGH (as processus for verification are mature) Cost of correction -CF1M1: LOW. The correction is basically a modification of a constant in the PWM management code. -CF1M2: LOW (same as CF1M1) -CF2M1: HIGH to VERY HIGH For the "silence time underestimation" fault model, the combination of a VERY HIGH cost, a LOW detection probability and a LOW correction cost leads the introduction of a new detection means. This means shall have a MODERATE cost and HIGH detection coverage.

  no short-cut Dissipated power < P max Application High / low ratio of PWM signal over Δ𝑇 < MR max mean speed over Δ𝑇 < MS max Contract is propagated Contract is propagated again Contract is undertaken

Figure 4 :

 4 Figure 4: Motor drivers contracts Another contract defined during the PWM controller design binds the hardware and software. It concerns the static and dynamic definition of the interface (see property C4). This contract is verified at the LT and AT modelling levels. Section 5.b describes the use of the contract verification on the PWM controller.

Figure 5 :

 5 Figure 5: Tool organization

Figure 7 :

 7 Figure 7: Functional model of Motor ControlThe PWM controller is integrated in the driver abstraction SystemC module by implementing two specific C++ methods for the definition of the API (init_mc_driver() and put_mc_pwm(int32 value)). The SystemC module integrates the driver abstraction with the PWM controller and the motor regulation. It declares two public methods init_mc() and put_setpoint(value) to allow access to PWM controller and to write into the setpoint data. It includes one thread (or method) activated every 10 ms calling c_mc_ctrl() for regulation control activation. Using the toolset, the SystemC motor control model is first wrapped up into a python module which is then imported into the tool environment as a new component of the virtual platform. To use this newly defined component, the user first instantiates and initializes it. Then, he/she sets the setpoint value using to the python API (driver_obj=vlab.get_instance("MC").obj, driver_obj.init_mcr(), driver_obj.put_setpoint(value)) .

Figure 8 :

 8 Figure 8: Test bench environment To observe the results of the P5 verification, a trace is placed on the P5 checker diagnosis output port using tool tracing capabilities. See the trace of the scenario result in Figure 11. Loosely Timed model The Loosely Timed Model is a structural and behavioral refinement of the Functional Model. It introduces several hardware and software components: -The timer IP abstraction generates the clocks used to generate the PWM signal. This IP is modelled in systemC. It integrates all design constraint defined in C4 property (register size, frequency range). -The driver abstraction provides the interface to the timer IP. This software component is defined in C. It is interfaced with the timer abstraction IP via a Hardware Software Interface (HSI). It integrates the C4 HSI contract (static and sequence interface). -A SystemC bus driver interface allows to access to the timer IP peripheral. The C to SystemC interface complies with the TLM2 loosely timed standard. It is available in the tool component library.

Figure 9 :

 9 Figure 9: LT model of the Motor Control

Figure 10 :

 10 Figure 10: AT Model of the Motor ControlThe test bench is only slightly adapted with respect to the one used at the more abstract levels. A python transactor component is inserted and connected to the same test bench driver. It integrates one method for target memory introspection using debugging API feature of the tool environment (write_memory()) and a second one to dummy the motor control initialization as it is now performed by the µC start-up code. The simulation scenario generator and the property checker component do not change. The P5 checker can be reused for the validation of the contract defined in the FMECA process. The functional property P4 of the PWM control is verified. The C4 contract on the HSI register access and sequence is also verified. The simulation trace typical of every refinement level is depicted below in the Figure11.

Figure 11 :

 11 Figure 11: Trace of dead zone validation

Seconded from Continental Automotive France, 1 avenue Paul Ourliac, 31036 Toulouse France

Seconded from Thales Avionics, 105 avenue du Général Eisenhower, 31100 Toulouse, France

Seconded from ACTIA Automotive, 5 rue Jorge Semprun, 31400 Toulouse, France

Acknowledgement

The authors thank all people and industrial partners involved in the Ingequip project. This work is supported by the French Research Agency (ANR) and by the industrial partners of IRT Saint-Exupéry Scientific Cooperation Foundation (FCS): Actia Automotive, Airbus, Airbus Defence and Space, ASTC Design Partners, Continental Automotive, SAGEM, Systerel, Space Codesign Systems and Thales Avionics.