Anthony Robert
email: anthony.robert@ifpen.fr

K Truffin

N Iafrate

S Jay

O Colin

C Angelberger

Large-eddy simulation analysis of knock in a direct injection spark ignition engine

Keywords: Knock, LES, Direct injection, Downsized engine, MAPO

Downsized spark ignition (SI) engines running under high loads have become more and more attractive for car manufacturers because of their increased thermal efficiency and lower CO2 emissions. However, the occurrence of abnormal combustions promoted by the thermodynamic conditions encountered in such engines limits their practical operating range, especially in high efficiency and low fuel consumption regions. One of the main abnormal combustion is knock, which corresponds to an autoignition of end gases during the flame propagation initiated by the spark plug. Knock generates pressure waves which can have long term damages on the engine, that is why the aim for car manufacturers is to better understand and predict knock appearance. However an experimental study of such recurrent but non-cyclic phenomena is very complex, and these difficulties motivate the use of CFD for better understanding them.

In the present paper, Large-Eddy Simulation (LES) is used as it is able to represent the instantaneous engine behavior and thus to quantitatively capture cyclic variability and knock. The proposed study focuses on the LES analysis of knock for a direct injection SI engine. A spark timing sweep available in the experimental database is simulated, and 15 LES cycles were performed for each spark timing. Wall temperatures, which are a first order parameter for knock prediction, are obtained using a conjugate heat transfer study. Present work points out that LES is able to describe the in-cylinder pressure envelope whatever the spark timing, even if the sample of LES cycles is limited compared to the 500 cycles recorded in the engine test bench. The influence of direct injection and equivalence ratio stratifications on combustion is also analyzed. Finally, focusing on knock, a MAPO (Maximum Amplitude Pressure Oscillation) analysis is conducted for both experimental and numerical pressure traces pointing out that LES well reproduces experimental knock tendencies.

Introduction

In the last decade, downsized spark ignition (SI) engines running under high loads behave attracted increasing interest thanks to their increased thermal efficiency and low CO2 emissions. However, the related high engine loads induce severe thermodynamic conditions in the combustion chamber, promoting the occurrence of abnormal combustions like knock. This phenomenon is related to an uncontrolled auto-ignition (AI) of the fresh gases before their consumption by the premixed flame initiated by the spark ignition device. It highly depends on the current combustion velocity as well as species composition and temperature fluctuations in the cylinder, making knock a recurrent but non-cyclic phenomenon. The detailed study of knock is experimentally complex because of the confined geometry and their possible destructive characteristics. Such limitations motivate the use of CFD for better understanding these specific combustions. Notably, RANS simulations were used to predict and understand the occurrence of knock (1; 2). However, the RANS approach is limited to the description of the mean cycle, which is not always affected by knock because of its sporadic nature. Large-Eddy Simulation (LES) appears as an attractive alternative because it allows the simulation of individual cycles, thus reproducing the link between knock and cyclic variability. In the last few years, first attempts were made to model these phenomena (3; 4; 5), demonstrating the potential of LES to address this topic. More recent work addressed first quantitative comparisons of LES results with experimental data [START_REF] Robert | LES prediction and analysis of knocking combustion in a spark ignition engine[END_REF], pointing out that LES is able to predict knock tendencies. LES also allows analyzing in detail the occurrence of knocking events. Low knock intensity is linked to one or several local auto-ignition spots which consume only the surrounding fresh gases, whereas the occurrence of a transition from deflagration to detonation is responsible for the highest knock intensities observed in SI engines [START_REF] Robert | LES study of deflagration to detonation mechanisms in a downsized spark ignition engine[END_REF]. The present paper addresses for the first time a quantitative LES study of knocking combustion in a SI engine, accounting for the direct injection and accurate wall temperature boundary conditions. A spark timing variation available in the experimental database is simulated in LES. The effect of injection on the mixture formation and in particular on equivalence ratio heterogeneity at spark timing is studied. The flow around the spark plug is then analyzed and observations are linked to the combustion behavior. Finally, the knock tendency as a function of the spark timing predicted by LES is compared with experimental trends.

EXPERIMENTAL SET-UP

The original engine configuration is a three cylinder fourvalve engine from groupe PSA, and the computational domain corresponds to one single cylinder chosen amongst them, as shown in Fig. 1.

The configuration is a direct injection spark ignition engine, whose features are summarized in Tab. 1. The operating point is at 5500 rpm with an IMEP of 19 bars. The fuel is a commercial European gasoline (SP95-E10) with 10% of volume fraction of ethanol. This fuel is injected thanks to an asymmetric five holes injector using a single injection occurring early during the intake phase in order to obtain a mixture as premixed as possible at top dead center (TDC). A spark timing (ST) sweep is available in the experimental database, and for confidentiality issues, spark timing values are given in this article relative to the latest spark timing called "reference spark timing".

MESHES AND BOUNDARY CONDITIONS FOR LES

Seventy tetrahedral meshes are needed to simulate a full cycle, owning between 1.7 and 9.2 Million cells. The LES methodology used for the moving mesh is described in [START_REF] Truffin | Using largeeddy simulation and multivariate analysis to understand the sources of combustion cyclic variability in a sparkignition engine[END_REF]. The mean cell size is approximately 0.05mm at the residual lift for intake and exhaust valves. Cell size in the chamber is around 0.7mm during the main part of the cycle and 0.5mm during combustion phases. Finally, the mesh size is refined down to 0.2mm in the vicinity of the spark plug during ignition. Several meshes are presented in Fig. 2.

At the inlet and outlet of the LES domain, boundary conditions are handled using the NSCBC approach (9; 10). The same temporal signals resulting from a GT-power calculation are imposed as boundary conditions for each LES cycle.

Wall temperatures, which are a first order parameter to study knock, are estimated using a conjugate heat transfer (CHT) based on RANS calculations. Computed temperature fields are imposed as wall boundary conditions in LES. Fig. 3 presents the resulting normalized wall temperature distribution for the cylinder head and the intake pipe. The hot region visible on the cylinder head is thus taken into account. Owing to the strong thermal inertia, these wall temperatures were kept constant for all simulated LES cycles. The valves were not included in CHT study, and their temperatures are estimated using thermal images of the real engine test bench (Fig. 4). They are then imposed in the LES calculations with a temperature gradient on the edges of the valves (Fig. 5) to approximate the heat exchange at the valve seats during valves closure.

NUMERICAL SET-UP

Large-Eddy Simulations are performed using the AVBP code (8; 11), co-developed and co-owned by CERFACS and IFPEN. AVBP solves the multi-species, compressible, reactive Navier-Stokes equations on unstructured and moving meshes [START_REF] Moureau | Numerical methods for unsteady compressible multi-component reacting flows on fixed and moving grids[END_REF]. Time advancement is explicit and convection is discretized in the present simulations using a second-order (in space and time) centered finite volume Lax-Wendroff scheme [START_REF] Lax | Systems of Conservation Laws[END_REF]. The subgrid scale turbulence is modeled by a Smagorinsky model [START_REF] Smagorinsky | General circulation experiments with the primitive equations. I: The basic experiment[END_REF] with a constant set to Cs=0.18. Finally, the wall shear stress and heat flux are imposed using a RANS modeling based on a logarithmic law-of-the-wall [START_REF] Schmitt | Large-eddy simulation and experimental study of heat transfer, nitric oxide emissions and combustion instability in a swirled turbulent high-pressure burner[END_REF].

The Lagrangian formalism [START_REF] Garcia | Comparison between Euler/Euler and Euler/Lagrange LES approaches for confined bluff-body gas-solid flow prediction[END_REF] is used to represent the liquid phase. Due to the diluted liquid phase assumption, droplets injected are smaller than flow characteristics scales and the point source approximation is employed. To initiate the spray physics at the nozzle exit, a simple disk model is used [START_REF] Iafrate | A study of gasoline multi-hole injector spray dynamics using LES[END_REF]. This model assumes that injected liquid can be described with a discrete approach omitting the dense region of the spray. A full disk surface is defined at the nozzle exit, on which a random deposit of particle is made. The main idea is to define a Gaussian velocity profile in order to conserve momentum. The following equation is verified:

V mean is computed from the instantaneous injection rate

Q inj : V mean (t) = Q inj S inj α l ρ p (3)
Where α l is the cavitation coefficient. The Gaussian velocity profile V (r) is defined by:

V (r) = V max exp(-σr 2 ) (4)
The maximal velocity is located at the center of the disk, and is:

V max = σV mean 1 -exp(-σ) (5)
To respect physical properties, the maximum velocity is limited by the Bernoulli velocity. As a consequence, the Gaussian parameter σ is adjusted to conserve momentum (Fig. 6).

Finally a turbulent fluctuation at the nozzle exit is introduced by adding a random velocity contribution to the normal and radial liquid velocity, as:

V (r, θ) = V (r)(1 + 2V rms * RN ) ( 6
)
Where θ is the angular position on the injection disk surface, V rms a given dimensionless velocity fluctuation and RN a random number between 0 and 1 different for each position (r,θ).

The flame front propagation is described using the ECFM-LES premixed combustion model (8; 17; 18). The local mean laminar flame speed is obtained from the correlation of Yahyaoui et al. [START_REF] Yahyaoui | IDETHANOL project : corrélation des vitesses de flamme laminaires des mélanges essence/éthanol[END_REF] that was developed for various mixtures of gasoline/ethanol. The spark ignition is modeled using the ISSIM-LES model [START_REF] Truffin | A spark ignition model for large eddy simulation based on an FSD transport equation (ISSIM-LES)[END_REF] which includes the description of the electrical circuit and 3D modeling of the flame kernel growth during ignition. The tabulated auto-ignition model (TKI-LES) (6; 21) is used to model the reaction rate linked to the occurrence of auto-ignition in the fresh gases. This model is based on a progress variable that is independent of the premixed flame one. Both phenomena are therefore fully decoupled so that an auto-ignition spot does not create artificially a propagating flame. The TKI model uses a look-up table built from complex chemistry simulations using a TRF surrogate with 42.8% isooctane, 13.7% nheptane, 43.5% toluene and the LLNL kinetic mechanism considering 1388 species and 5935 reactions [START_REF] Mehl | Kinetic modeling of gasoline surrogate components and mixtures under engine conditions[END_REF]. Based on local conditions, an auto-ignition reaction rate is extracted from the TKI table for each node of the mesh, and used in the species transport equations.

The post-flame kinetics is taken into account to correct the burned gases state and temperature. First the remaining fuel that has not been consumed by the propagating flame can be post-oxidized using the following consumption rate model:

ωF b = ρW F A ′ exp - E ′ a RT b ρ Y b F W F 0.55 ρ Y b O2 W O2 0.9 (7
) with A ′ = 6.10 11 cm 3 .mol -1 .s -1 and E ′ a = 41500 cal.mol -1 . The b exponent stands for the burnt gas state.

Second the kinetic oxidation of CO is also introduced:

CO + 1 2 O 2 ↔ CO 2 (8)
The rate of reaction of Eq. 8 is given by:

Q = k f ρ Y b CO W CO ρ Y b O2 W O2 1/2 -k r ρ Y b CO2 W CO2 (9)
with W CO , W CO2 and W O2 the molar mass of CO, CO2 and O2 respectively, and

Y b CO , Y b CO2 , Y b
O2 their mass fractions in the burnt gases. k f and k r are the forward and reverse reaction rates:

k f = A exp - E a RT b and k r = k f K eq ( 10
)
where the pre-exponential constant A is fixed to 9.10 8 cm 3 .mol -1 .s -1 , the activation energy is 40000 cal.mol -1 and K eq is the equilibrium constant defined by [START_REF] Kuo | Principles of combustion[END_REF]. The rate constants of this reduced two-step mechanism were adjusted to recover the correct heat release rate and burned gases temperature on homogeneous auto-ignition reactor calculations at various fuel/air ratios and pressures.

For this purpose the same LLNL mechanism used to generate the TKI table was used as a reference solution.

RESULTS AND DISCUSSIONS

Fifteen consecutive cycles for the reference spark timing (the latest one of the ST sweep) were first simulated with LES.

Reference spark timing analysis

The simulated operating point owns an early injection illustrated on Fig. 7. Fuel is injected during intake valve opening, droplets are driven into the tumble motion and this type of injection should lead to a premixed mixture at TDC. To analyze injection in more details, Fig. 8 shows the temporal evolution of the liquid mass in the combustion chamber, which increases during injection before reaching a maximum value around -200 CAD before TDC. Liquid masse then decreases due to the evaporation of droplets which are fully evaporated around 60 CAD bTDC. Fig. 9 shows the distribution of local fuel/air equivalence ratio in the chamber at three different instants for the simulated cycles. The mixing resulting from the compression of the tumbling intake flow, reduces the spatial variations as one comes close to TDC, where it centers around the targeted mean equivalence ratio of 1.2.

However, the mixture is not fully homogeneous at TDC, exhibiting a non-negligible instantaneous stratification, and shows important cycle-to-cycle variability. Fig. 10 shows the spatial distributions for cycles 2 and 10 at TDC, which are two extreme ones in terms of mixture heterogeneity in the LES cycles. The spatial variance of equivalence ratio for cycle 10 is higher than for cycle 2, with equivalence ratio extrema going from 0.9 to 1.6, which is much more heterogeneous than cycle 2. It is important to notice that cycle 2 is one of the most homogeneous LES cycles. These observations are confirmed on Fig. 11 where equivalence ratio on a horizontal cut plane at TDC is presented for the same LES cycles. The spatial distribution of equivalence ratio is much smaller for cycle 2 than for cycle 10. Overall, the LES predictions indicate that a perfectly premixed fuel/air mixture is not obtained at TDC.

Despite the found strong differences in mixture homogeneity between cycles 2 and 10, this is not sufficient to explain the predicted cyclic combustion variability. As shown in Fig. 12, these two extreme cycles in terms of heterogeneity are indeed not extreme in terms of combustion speed, as both are located in the bottom part of the cylinder pressure envelope. Overall the cyclic combustion variability predicted by the 15 LES cycles reproduces quite well the one observed experimentally based on 500 cycles. Note that the reported pressure was recorded at the same position in experiments and LES.

However, it may be noticed that LES cycle number 15 lies below the statistically probable envelope of the experiments, right after the spark ignition timing. To understand what happens in more details, Fig. 13 and Fig. 14 show the mean equivalence ratio and flow velocity magnitude in a sphere of 5 mm around the spark plug at the instant of spark ignition. Values are plotted over the cycle identification number, and we can first notice that cycle 1 is not taken into account as results depend too much on the initial condition.

The mean equivalence ratio is fluctuating around the mean value of 1.2, with limited cycle-to-cycle variability. In particular the mean equivalence ratio at the spark of cycle 15 is very close to the one for cycle 11 which is the second fastest cycle predicted in LES (see Fig. 12). On the other hand, cycle 2 and 10 have different mean equivalence ratios at the spark plug but exhibit very similar combustion speeds. It appears that despite the fact that equivalence ratio at the spark plug is one of the parameters influencing the ignition phase and combustion development, its cyclic variability does not appear to be sufficient to explain the combustion variability. Fig. 14 shows that the cycle-to-cycle variability of mean flow velocity at the spark is important, with a factor of 2 between the maximum and minimum level. The cycle exhibiting the smallest mean velocity is cycle 15, which is the slowest combustion cycle. Cycles 2 and 10 exhibit very similar levels, but lowers than the one predicted for cycle 11, which is one of the fastest cycles. Fig. 15 compares isosurfaces of premixed flame progress variable colored by temperature for cycles 11 and 15, at 16 CAD after spark ignition. A vortex pushes the flame kernel of cycle 15 down towards the piston. As a result, the ensuing propagation appears to be slowed down as compared to cycle 11. In the latter, a high horizontal flow velocity favors the horizontal propagation of the flame, which interacts much later with the piston. All these observations points out that ignition and flame propagation are more influenced by the velocity magnitude at the spark plug and the direction of the flow in the combustion chamber than by cyclic variability of equivalence ratio.

Methodology for simulating the spark timing sweep

A spark timing sweep covering eight ignition times is available in the experimental database. Six of them are simulated in LES, respectively -7.51 CAD, -6.76 CAD, -5.26 CAD, -3.76 CAD, -2.25 CAD before the reference ST and the reference one. The strategy followed to perform the LES study of the ST sweep is illustrated in Fig. 16. In an approach already proven accurate in the past [START_REF] Robert | LES prediction and analysis of knocking combustion in a spark ignition engine[END_REF], the fifteen full consecutive LES cycles at the reference ST are used to yield initial conditions just before ignition for the studied ST. As a result only the combustion phases need to be computed. Several cycles can thus be simulated in parallel allowing reducing return times.

This approach is justified in the absence of dependency between consecutive cycles, which was verified for the reference spark timing of the studied engine and supposed to be valid also for all studied ST. This assumption may not be valid in the general case.

Analysis of the spark timing sweep

Fig. 17 to Fig. 19 compare the in-cylinder pressure predicted by the 15 LES cycles (black lines) with the 500 experimental cycles (brown lines) for three of the simulated ST. The LES cycles are able to reproduce the experimental pressure envelope recorded at the cylinder head whatever the spark timing (results for the ST not presented here have the same good agreement).

To go further on the analysis of in-cylinder pressure, the evolution of the coefficient of variation (COV) of the maximum pressure is presented on Fig. 20 for the whole spark timing sweep. The COV is defined as the ratio of the standard deviation to the mean of the maximal in-cylinder latest spark timing is well reproduced by the LES, even if the numerical results overestimate a bit experimental ones, mainly due to the limited sample of LES cycles.

To compare combustion velocities between LES and experiment, the time ∆t to go from 10% (CA 10 ) to 90% of burned mass in the combustion chamber (CA 90 ) is computed:

∆t = CA 90 -CA 10 (11)
Fig. 21 shows the resulting mean error ǫ between LES results and experimental findings for the simulated ST, computed as:

ǫ = ∆t Exp.mean -∆t LESmean ∆t Exp.mean (12)
where the respective mean values are computed as an ensemble average of the available cycles (15 in LES, 500 in experiments). The maximum error is around 8%, confirming a good reproduction by LES of the experimentally observed combustion speeds. LES is found to underestimate it for the early ST, and to overestimate it for the latest spark timing. The reasons for this have not yet been further investigated.

Focusing on knock, the use of models ECFM-LES and TKI-LES allows to follow distinctly the propagation of the premixed flame initiated at the spark plug and the autoignition spots. Fig. 22 illustrates the knock occurrence for the LES cycle n o 12 at a RST of -7.51 CAD. Premixed flame position is computed using a contour of ECFM progress variable colored by temperature, whereas autoignition spot is represented by a green contour of TKI progress variable. For this cycle, the premixed flame propagation is slow down in a rich region (equivalence ratio field not shown here) under exhaust valves giving much more time to fresh gases to autoignite. Such analysis allows to validate the occurrence of autoignition in the engine but can only be done for numerical results. As the aim here is to compared LES knock results to experimental ones, and because knock is characterized by pressure oscillations, fluctuation of pressure traces available both for LES and experiment are now analyzed.

From Fig. 12 and Fig. 17 to Fig. 21, the pressure fluctuation amplitudes increase when spark timing occurs more and more early during the cycle. To quantify this qualitatively accurate tendency predicted by LES, a Maximum Amplitude of Pressure Oscillations analysis currently (MAPO) [START_REF] Zhen | The engine knock analysis -An overview[END_REF] is conducted. This analysis is based on the pressure signal, the same numerical treatment being applied to LES and experimental signals. First, a sensitivity study of the MAPO results on the location of the pressure sensor is conducted using LES signals. To this purpose, three other numerical sensors are located between the valves as presented in Fig. 23.The MAPO analysis is conducted on these four probes for the fifteen LES cycles at the reference spark timing. Fig. 24 shows the resulting dependency of the knock intensity on the 15 cycles as a function of the position pressure recording. It reveals the existence of two groups with different values of knock intensity. Sensors n o 1 and n o 3 show equivalent knock intensity, much higher than knock intensity level predicted with sensors n o 2 and n o 4. This observation underlines the important sensitivity of pressure sensor location on the cylinder head.

In what follows, the MAPO analysis concerns the signals recorded at sensor n o 3, corresponding to the experimental location. In Fig. 25 to Fig. 27, the knock intensity is plotted over the cycle identification number. To make these plots more readable, and as the LES sample is much smaller than for experiments, the identification number for LES cycles is multiplied by 10 in these figures: so LES cycle 2 becomes cycle 20, LES cycle 3 become cycle 30 ... etc. Mean values are also plotted both for LES and experiment. For all spark timings, the knock intensity predicted by LES is in good agreement with the experimentally observed levels. Based on the mean knock level evolution, the earlier spark ignition occurs in the cycle, the higher the knock intensity is, both in the experiments and in the LES. Results are also consistent for the three other spark timings simulated in LES, but not presented in this paper. Another way to quantify the knock tendencies is to look at the percentage of knocking cycles. A cycle is considered as knocking when its knock intensity is superior to a threshold value defined on the real engine test bench. For the operating conditions of this study, the value of 1.1 bar is used both for LES and experiment. Fig. 28 shows the evolution of the percentage of knocking cycles over the relative spark timing (RST). As the sample of LES cycles is reduced, the experimental sample is divided randomly into subsets of 15 cycles for which the analysis is repeated, in addition to the full 500 cycles sample. The subset yielding the highest percentage is also plotted. At the reference ST (RST equal to zero), non knocking cycles are observed both for LES and experimental results. First knocking cycles are detected at -2.25 CAD before the RST in the experiment and at -3.76 CAD in the LES, with very comparable intensity level (around 5%). The main difference occurs at a RST of -5.26 CAD where the experimental level falls down close to zero, whereas LES predicts a continuous increase of the percentage of knocking cycles. The experimental behavior is quite puzzling and difficult to explain, but was not further explored here.

Two additional earlier ST have been simulated in LES, but weren't studied experimentally to avoid damages. The evolution of the percentage of knocking cycles is presented in Fig. 29. After a constant level until a RST of -11.3 CAD, LES predicts a sharp increase of the number of knocking cycles to reach a value close to 90%. This observation justifies the fact the real engine should not run with such early ST, for which the knock intensity (not shown here) increases sharply.

CONCLUSIONS

This paper presents for the first time a quantitative LES study of cyclic variability and knock taking into account direct fuel injection in a production four-valve spark ignition engine fueled with a typical European gasoline, and running at high load and speed. Particular attention has been paid to an accurate estimation of wall temperatures using a CHT study based on RANS calculations. Results of this CHT study were imposed as boundary conditions in LES. First, an analysis of the in-cylinder mixture preparation as a result of fuel injection and intake flow revealed that despite an early injection timing, the tumble motion is not able to mix air and fuel sufficiently fast to obtain a perfectly premixed mixture at spark timing. LES was shown to be able to reproduce the experimental cylinder pressure variation envelope using a limited sample of 15 simulated cycles, and this for all studied spark timings. Even if cycle-to-cycle variability is well predicted by LES, one LES cycle was found to lie below the experimentally reported envelope. A detailed analysis of the mean flow around the spark plug pointed out that cyclic equivalence ratio variability is limited and was probably not the main reason for cycle-to-cycle combustion variability. It was shown that ignition and flame propagation were more influenced by the velocity magnitude at the spark plug, and by the cyclic variability of the flow during early flame propagation. To quantify knock intensity, the same MAPO analysis was applied to LES and experimental results. A sensitivity to the pressure sensor location on the cylinder head was conducted and showed that it strongly influenced the resulting knock intensity. This confirms that this location should be chosen with care for an on-board knock detection. The dependency of knock on spark timing was then analyzed. LES was shown to reproduce quite accurately the experimental findings in terms of intensity and percentage of knocking cycle's variation with spark timing. Future work concerns the exploration of the causes of cyclic variability and knock, in an attempt to identify control parameters allowing to reduce them in early design phases.

Figure 1 .

 1 Figure 1. Presentation of the real engine configuration and the LES computation domain.

  Spark timing sweep From -7,5 CAD before the reference spark timing to the reference ST

Figure 2 .

 2 Figure 2. Presentation of meshes used during intake valve opening (left and top right) and at the spark plug during ignition (bottom right).

Figure 3 .

 3 Figure 3. Normalized wall temperatures imposed as LES boundary conditions for cylinder head and intake pipe.

Figure 4 .

 4 Figure 4. Thermal image of the exhaust valves.

Figure 5 .

 5 Figure 5. LES wall temperatures of exhaust valves.

Figure 6 .

 6 Figure 6. Example of adapted normalized velocity profiles (top) depending on the instantaneous injected mass flow rate (bottom)

Figure 7 .

 7 Figure 7. Image of injected droplets in the cylinder during intake valve opening for one LES cycle. Droplets are colored by their temperature.

Figure 8 .

 8 Figure 8. Temporal evolution of liquid mass for the fifteen LES cycles at the reference ST.

Figure 9 .

 9 Figure 9. Spatial distribution of equivalence ratio for fifteen LES cycles at three instants during compression.

Figure 10 .

 10 Figure 10. Spatial distribution at TDC for two individual cycles at the reference ST.

Figure 11 .

 11 Figure 11. Comparison of equivalence ratio fields on a horizontal cut plane at TDC for two LES cycles.

Figure 12 .

 12 Figure 12. Normalized in-cylinder pressure comparison between LES (black) and experiment (brown) at the reference spark timing.

Figure 13 .

 13 Figure 13. Mean equivalence ratio at the spark plug for the fifteen LES cycles at the reference spark timing.

Figure 14 .

 14 Figure 14. Mean velocity magnitude at the spark plug for the fifteen LES cycles at the reference spark timing.

Figure 15 .

 15 Figure 15. Evolution of the flame kernel at 16 CAD after spark ignition for the slowest (cycle 15) and the fast (cycle 11). Arrows indicates the main direction of flame propagation.

Figure 16 .

 16 Figure 16. Illustration of the methodology used to simulate the spark timing sweep. Only combustion phases are simulated for the different spark timings, starting from the aerodynamic fields obtained from the reference case.

Figure 17 .

 17 Figure 17. Temporal evolution of normalized in-cylinder pressure at a relative spark timing (RST) of -2.25 CAD.

Figure 18 .

 18 Figure 18. Temporal evolution of normalized in-cylinder pressure at a RST of -5.26 CAD.

Figure 19 .

 19 Figure 19. Temporal evolution of normalized in-cylinder pressure at a RST of -7.51 CAD.

Figure 20 .

 20 Figure 20. Evolution of the coefficient of variation (COV) of the maximum pressure on the whole spark timing sweep.

Figure 21 .

 21 Figure 21. Combustion duration error over the relative spark timings.

Figure 22 .

 22 Figure 22. Autoignition spot (green) for LES cycle n o 12 at +18 CAD after TDC. Premixed flame position is computed using a contour of ECFM progress variable colored by temperature.

Figure 23 .

 23 Figure 23. Localization of the four numerical sensors on the cylinder head.

Figure 24 .

 24 Figure 24. MAPO analysis of the pressure signal recorded at different sensor locations.

Figure 25 .

 25 Figure 25. Knock intensity comparisons between LES and experiment at the reference ST. Mean knock intensity level is 0.34 bar for experiment and 0.37 bar for LES.

Figure 26 .

 26 Figure 26. Knock intensity comparisons, ST is set to the reference minus 3.76 CAD. Mean knock intensity level is 0.5 bar for experiment and 0.59 bar for LES.

Figure 27 .

 27 Figure 27. Knock intensity comparisons, ST is set to the reference minus 7.51 CAD. Mean knock intensity level is 0.64 bar for experiment and 0.85 bar for LES.

Figure 28 .

 28 Figure 28. Percentage of knocking cycles over the relative spark timings for the common ST between experiment and LES.

Figure 29 .

 29 Figure 29. Percentage of knocking cycles over the relative spark timings for all the ST simulated in LES.

Table 1 .

 1 Engine characteristics and simulated operating point.

Prepared using sagej.cls

Acknowledgements

This work was funded by the GSM (Groupement Scientifique Moteur: groupe PSA, Renault SAS and IFP Energies nouvelles) and by ADEME (French Agency for Environment and Energy Control). The authors want to acknowledge Benoit ENAUX from groupe PSA for providing engine database and CHT calculation. NOMENCLATURE AI : auto-ignition TDC: Top Dead Center ST : Spark timing RST : Relative Spark Timing CHT : Conjugated Heat Transfer MAPO : Maximum Amplitude of Pressure Oscillations


Marc Vincent

Olivier Mathieu

Patrick Nolain

Cecilia Menacé

Dr Sonia Khier
email: sonia.khier@umontpellier.fr

Population pharmacokinetics of levobupivacaine during a transversus abdominis plane block in children

come

Introduction

Over the past years transversus abdominis plane (TAP) block, a loco regional anaesthesia technique, has become increasingly popular in abdominal surgery because it provides very effective pain relief both pre-and post-operatively thereby reducing postoperative opioid requirement. A local anaesthetic solution is injected into a neurovascular space, between the transverse and internal oblique muscle of the abdomen. The solution anesthetizes the nerves supplying the anterior abdominal wall (T9 to L1) [START_REF] Young | Clinical implications of the transversus abdominis plane block in adults[END_REF][START_REF] Willschke | Pediatric regional anesthesia: abdominal wall blocks[END_REF], Figure X. Levobupivacaine is commonly used as an anaesthetic in TAP block. The solution diffuses into the nerve, and causes a reversible block of the sodium channels. This prevents depolarisation and conduction of the nerve action potential. Levobupivacaine local concentration is correlated to its efficiency. But the solution is also absorbed in the systemic and may induce adverse effects (carditoxicity, neurotoxicity) [START_REF] Mcleod | Levobupivacaine[END_REF]. Absorption from the site of administration is affected by the vascularity of the tissue. Protein binding is upper to 97%. Levobupivacaine is extensively metabolized with no unchanged levobupivacaine detected in urine or faeces. CYP3A4 and CYP1A2 mediate the metabolism of levobupivacaine to desbutyl-levobupivacaine and 3-hydroxy levobupivacaine, respectively The dosing regimen is still empirical in paediatric surgery and rarely based on pharmacokinetic data. In practice, an injected volume of 0.2 à 0.3 ml/kg of levobupivacaine (2.5 mg/ml) is used [START_REF] Fredrickson | Improved analgesia with the ilioinguinal block compared to the transversus abdominis plane block after pediatric inguinal surgery: a prospective randomized trial[END_REF][START_REF] Bhalla | Ultrasound-guided trunk and core blocks in infants and children[END_REF].

The aim of this work was therefore to describe pharmacokinetics of levobupivacaine and investigate possible sources of between subject variability in children by a population approach when performing an ultrasound-guided TAP block in ambulatory abdominal wall surgery. As the volume of anaesthetic injected could possibly influence its pharmacokinetics [START_REF] Hong | A comparison of high volume/low concentration and low volume/high concentration ropivacaine in caudal analgesia for pediatric orchiopexy[END_REF], the study was conducted according to two methods of administration: high or low volume of local anaesthetic injected.

Methods

This prospective, randomized, double-blind study protocol was approved by the institutional ethical committee (CPP Sud Mediterranean IV, 2013-A00618-37, Montpellier, France), and registered in Clinical Trials Database (ClinicalTrials.gov. NCT02064088). The written informed consent was obtained from the parents or legal guardians of all participating subjects before randomization.

Patients

Children from 1 to 5 years old (14 to 71 months) and planned for inguinal surgery were candidate to the study. Exclusion criteria included any contraindication to general anaesthesia or TAP block, amide local anaesthetic drug allergy, American Society of Anesthesiologists physical status 3 and 4, nonsteroidal anti-inflammatory drug intolerance, known cardiac / renal / hepatic dysfunction, parental refusal, and concurrent participation to another clinical trial.

Study design and sampling protocol

Each patient received 0.4 mg/kg of levobupivacaine (Chirocaine®, Abbott, Rungis, France) [START_REF] Sola | Ultrasound-guided tranversus abdominis plane block for herniorrhaphy in children: what is the optimal dose of levobupivacaine?[END_REF] and were assigned to one of the two study groups using a computer generated random sequence. In group 1 (20 patients), the volume of the local anaesthesia solution was adjusted to 0.2 mL/kg of 0.2% levobupivacaine (low volume/high concentration). In group 2 (20 patients), the volume of the local anaesthesia solution was adjusted to 0.4 mL/kg of 0.1% levobupivacaine (high volume/low concentration). Blood samples were collected at 5, 15, 20, 25, 30, 45, 60 and 75 min following the block injection. To quantify the unbound fraction of levobupivacaine, orosomucoid (α1-acid glycoprotein) plasma was immediately ultra-filtrated (cut-off 30 kDa, Ultracel®, Merck Millipore®, Tullagreen, Ireland) at room temperature and stored at -20°C until analysis. For plasma, the assay was done by simple protein precipitation in acetonitrile with added internal standard (etidocaine) and dilution in the mobile phase. The analysis of the prepared samples was performed by liquid chromatography combined with a hybrid triple quadrupole trap mass spectrometer (3200 Qtrap, ABSciex®, Les Ulis, France). Free levobupivacaine concentrations were analyzed by direct injection of the ultra-filtrated serum. According to a previously published method regarding ropivacaïne, the interand intra-day percentage coefficients of variation were below 8 % and 6%, respectively. The lowest limit of quantification (LLOQs) was 1 ng/ml [START_REF] Mathieu | Liquid chromatography-electrospray mass spectrometry determination of free and total concentrations of ropivacaine in human plasma[END_REF].

Pharmacokinetic analysis

Pharmacokinetic analysis was performed with total concentrations of levobupivacaine, by non-linear mixed-effects modelling using NONMEM® version 7.4 (ICON Development Solutions, MD, USA).

Based on a visual inspection of the individual pharmacokinetic profiles and a review of the literature, one-and two-compartment models were considered to describe the concentration-time data. Different kinetics of absorption (zero-or first-order) were investigated. The between-subject variability (BSV) for each model parameter was modelled using Equation 1:

𝑃𝑖 = 𝑃𝑝𝑜𝑝. 𝑒 𝜂𝑖 , 𝜂~𝑁(0, 𝜔 2 ) (1)
where Ρi is the value of the parameter "(e.g., ka or absorption rate constant, CL or clearance, V or volume of distribution) for the ith patient, Ppop is the value of the population parameter and η is the random variable. Usual error models (additive, proportional and combined) were evaluated to describe residual unexplained variability during the screening step. The BSV is modelled in terms of random effect (η) variables. Each of these variables is assumed to have mean 0 and a variance denoted by ω 2 , which is estimated. Population parameters (typical values and BSV) were estimated using the FOCE-I algorithm (First Order Conditional Estimates with Interaction).

The model selection was performed on the basis of the minimization status (i.e., successful) and the examination of standard error of the parameter estimates (i.e., < 50 %). Models were compared with regards to the objective function value (OFV = -2. log-likelihood). The best model was selected based on the lowest OFV and the inspection of goodness of fits plots, i.e., observed plasma concentrations versus population predicted concentrations (PRED) plots, observed plasma concentrations versus individual predicted concentrations (IPRED) plots and conditional weighted residuals (CWRES) versus time, PRED or IPRED scatterplots. Demographic characteristics (age, weight and gender) as well as the orosomucoide rate (Table 1) and the volume of the local anaesthesia solution injected (0.2 or 0.4 mL/kg) were evaluated as potential covariates. All potential covariates were added separately to the model (forward inclusion-backward elimination method) using different functions (linear, power, and exponential). All the functions were tested at each step of the building model. The significance of a parameter-covariate relationship was reflected by a decrease of OFV of at least 3.84 (p-value < 0.05, from a chi-squared distribution, with one degree of freedom). a: first quartile b: third quartile.

N

Model evaluation

The robustness of the model and the accuracy of parameters estimates were assessed using a bootstrap method [START_REF] Parke | A procedure for generating bootstrap samples for the validation of nonlinear mixed-effects population models[END_REF]. The entire procedure was undertaken in an automated fashion. From the original dataset, 1 000 bootstrap sets were drawn with re-sampling. For each of the 1 000 bootstrap sets, the pharmacokinetic population parameters were estimated. The corresponding mean, median, standard deviation, 2.5 and 97.5th percentiles were calculated. Visual predictive checks (VPC) [START_REF] Yano | Evaluating pharmacokinetic/pharmacodynamic models using the posterior predictive check[END_REF] were also performed to ensure that simulations based on model could reproduce the observed data. Perlspeaks-NONMEM Tools were used for both VPC and bootstrap evaluation [11].

Computation of individual exposures

Based on Empirical Bayes Estimates (EBE), single dose exposure values were calculated: ) (𝐾𝑎 -𝐾𝑒) ⁄ .

TAP effectiveness evaluation

The cutaneous incision was allowed after a minimum of 15 minutes post TAP injection. Insufficient block (TAP block failure) was defined as an increase of either heart rate or mean arterial blood pressure ≥ 20% at the skin incision, compared to reference values prior block placement. In case of insufficient block, remifentanil infusion was used as backup.

Statistical analysis

Standard statistical tests (t-test, independent Wilcoxon or Fisher's Exact Test) were performed to compare the two groups "TAP Successful" and "TAP Failed". All the pharmacokinetic parameters and exposure values were compared between the two groups. Statistical analyses were performed using R version 3.5.1. The median maximum concentration was 0.315 mg/L at 17 minutes after the injection and the median total drug exposure across time was 41 mg/L. min. Mean protein binding estimate was 99 %. All the individual exposures values are presented in Table 1. Individual pharmacokinetic profiles were obtained and presented in Supplementary Data.

Discussion

Levobupivacaine 0.25% can be administered with doses ranging up to 1 mL/kg [START_REF] Fredrickson | Improved analgesia with the ilioinguinal block compared to the transversus abdominis plane block after pediatric inguinal surgery: a prospective randomized trial[END_REF][START_REF] Bhalla | Ultrasound-guided trunk and core blocks in infants and children[END_REF][START_REF] Sola | Ultrasound-guided tranversus abdominis plane block for herniorrhaphy in children: what is the optimal dose of levobupivacaine?[END_REF][START_REF] Chen | Analgesic efficacy of transversus abdominis plane block in neonates and early infants for colostomy and reversal of colostomy[END_REF][START_REF]Assessment of the paediatric needs anaesthesiology[END_REF][START_REF] Sahin | Comparison of 3 different regional block techniques in pediatric patients. A prospective randomized single-blinded study[END_REF][START_REF] Sahin | Ultrasound-guided transversus abdominis plane block in children: a randomised comparison with wound infiltration[END_REF].These recommendations are rarely based on pharmacokinetic data. Previous studies with ropivacaine showed an influence of the injected volume on the effectiveness of regional anaesthesia [START_REF] Hong | A comparison of high volume/low concentration and low volume/high concentration ropivacaine in caudal analgesia for pediatric orchiopexy[END_REF][START_REF] Trifa | The analgesic effects of ropivacaine in ilioinguinal-iliohypogastric nerve block in children--concentration or volume? Middle East[END_REF] while no study investigated the influence of the injected volume of levobupivacaine during a TAP block. The purpose of this work was to describe the pharmacokinetics of levobupivacaine in children and investigate possible sources of between subject variability including volumes of solution injected (small volume/high concentration or high volume/low concentration) during TAP block.

The absorption rate from TAP to systemic circulation is fast, with a mean absorption half-life = 6 min.

BSV on Ka has been estimated at 81 %. The forward inclusion-backward elimination method did not identify a covariate for Ka, and particularly the injected volume considered as categorical covariate (high/low injected volume). So, from our results, the injected volume does not influence Ka.

Bioavailability factor (F) for "high volume injected" group is comparable to "small volume injected" group (F=0.88, [0.59-1.17]). It is important to note that whatever the group, the injected volumes were low (1.8 to 8.4 mL), which can lead more difficult to highlight difference. Apparent clearance (CL/F) was estimated at 0.15 L/min. This value is higher than the clearance reported for younger children during epidural block [START_REF] Chalkiadis | Age and size are the major covariates for prediction of levobupivacaine clearance in children[END_REF]. This is coherent with pharmacokinetics of levobupivacaine, metabolized by CYP3A4 and CYP1A2 isoforms. As the maturation of enzyme producing organs increases with age, an increase of mature cytochromes involve a higher clearance. We found a BSV on clearance and distribution volume similar to that reported by Chalkiadis et al. (epidural regional block, with a BSV for V and CL of 48.5% and 35.2% respectively) [START_REF] Chalkiadis | Age and size are the major covariates for prediction of levobupivacaine clearance in children[END_REF]. Our results confirm that weight must be considered when determining the dose in children.

The dosing regimen of 0.4 mg/kg is safe, no signs of toxicity were observed during the study. But an optimization can be done in terms of efficacy: TAP block effectiveness (15 min post injection) was inadequate for 13 patients out of 40. We have tried to establish a pharmacokinetic/pharmacodynamic (PKPD) model but PK parameters and efficacy were not correlated (Emax, Hill or Logit function). This finding is not surprising given that is an a posteriori analysis, the study was not designed for a PKPD analysis. The only PD endpoint was the increase of either heart rate or arterial blood pressure at the skin incision time (approximatively 15 min post injection), compared to reference values prior block placement. If the heart rate or arterial blood pressure increased, the TAP block was considered as failed (DV = 0). If they were stable, the TAP block was considered as successful (DV = 1). A PKPD model could not be establish without several and accurate PD endpoints available. We observed a significant difference in response to anaesthesia with a higher rate of success for girls (91.7%) than for boys (57.1%). One possible explanation for this difference could be related to weight. On average, girls had slightly higher weight than boys (15.7 kg for girls versus 14.6 kg for boys) and therefore received a higher dose into the deposit site (dose = 0.4 mg/kg). From the pharmacokinetic profiles, the plasma concentrations observed (Cmax = 0.315 mg/L and Cu,max = 3.15 µg/L, respectively) are lower than concentrations considered to be well tolerated for adults (2-4 mg/L and 110-300 µg/L, respectively [START_REF] Frawley | Pharmacokinetics of levobupivacaine following infant spinal anesthesia[END_REF]), even for an outlier patient (Cmax = 1.29 mg/L) (Figure 3). Considering these results, an adjustment of dosage could be considered for children, in order to improve efficiency of analgesia. or for failed TAP block (pink line).

Conclusion

The optimal dosing regimen of levobupivacaine for an effective TAP block has not been reported in the literature and is still explored [START_REF] Abdallah | Transversus abdominis plane block: a systematic review[END_REF]. In paediatrics, even though recommendations exist, the optimal dose to be administered is not clearly established and is still empirical. At 0.4 mg/kg, none of the 1 to 5 years-old children presented toxicity signs and there is some clues that this dose could be increased, to obtain a better efficacy. Finally children's weight must be considered to evaluate their clearance and thus to anticipate any risk of toxicity.

Fig. x

 x Fig. x Intraoperative transversus abdominis plane (TAP) blockade. Under ultrasound guidance, a regional block catheter is passed into the transversus abdominis plane at the anterior-axillary line, midway between the costal margin and iliac crest. An oblique cross-section at the level of T11 is shown (right). EOM, external oblique muscle; IOM, internal oblique muscle; TAM, transversus abdominis muscle. From Jablonka et al.

  𝑒 -𝐾𝑒.𝑡𝑚𝑎𝑥 -𝑒 -𝐾𝑎.𝑡𝑚𝑎𝑥 ), iii) time to reach Cmax : 𝑡𝑚𝑎𝑥 = 𝑙𝑛 ( 𝐾𝑎 𝐾𝑒

Fig. 1

 1 Fig. 1 Commonly goodness of fit plots for the levobupivacaine model are used for evaluating model misspecification. Observed concentrations versus: (A) Population Predicted concentrations or PRED, (B) Individual Predicted concentrations or IPRED. Conditional weighted residuals (CWRES) versus: (C) time, (D) PRED or (E) IPRED.

Fig. 2

 2 Fig. 2 Visual predictive checks (VPC) plot generated from 1 000 simulations. Pink areas represent the 95 % confidence intervals (CI) of the 5 th and 95 th percentiles of the predictions. Upper dotted line: 95 th percentile of the observations, lower dotted line: 5 th percentile of the observations. Blue area represents the 95 % CI of the median predictions. Middle dotted line: median of predictions, middle solid line: median of observations. Black dots represent observed plasma concentrations.

Fig. 3 (

 3 Fig. 3 (A) Individual concentration profiles of successful TAP block, (B) Individual concentration profiles of failed TAP block, (C) Profile of average concentrations for successful TAP block (blue line)

Acknowledgments

Thank you to Gregory Nickson, for the linguistic revision. This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Results

Pharmacokinetic analysis

A total of 247 sample points from 40 subjects were included in the dataset. All the concentrations under the low limit of quantification (LLOQ = 0.05 mg/L) were considered as missing values (MDV=1).

The data were best described by a one-compartment model parameterized in terms of an absorption rate constant (ka, h-1) characterizing the first-order absorption process from the depot to the central compartment, an apparent distribution volume (V/F, L) and clearance (CL/F). The residual unexplained variability was described with a proportional error model (Equation 2):

where Ci is the concentration observed to the ith time and Cpred the concentration predicted by the model.

Weight significantly decreased the objective function and BSV (from 50 to 41%) of clearance and was integrated as covariate: Population pharmacokinetic parameters are presented in Table 2. Figure 1 clearance, V: volume of distribution of central compartment, ka first-order absorption rate constant.

Population pharmacokinetic parameters

Pharmacodynamics of levobupivacaine

Among the 40 patients, none of the 1 to 5 years-old children presented toxicity signs. An insufficient TAP block anaesthesia was obtained for 13 patients and a total anaesthesia was obtained for 27 patients. Pharmacokinetic parameters were not significantly different between the two groups:

effective TAP block versus failed TAP block. PK profiles are presented in Supplementary Data. But we observed that gender and effectiveness of TAP block were dependent variables. (Fisher's Exact Test, p-value = 0.03).


Yoann Blein

Arnaud Clère
email: arnaud.clere@minmaxmedical.com

Fabrice Bertrand
email: fabrice.bertrand@blue-ortho.com

Yves Ledru

Roland Groz

Lydie Du Bousquet

Improving Trace Generation and Analysis for Medical Devices

Introduction

A new generation of medical devices emerges to support increasingly more complex medical decisions and procedures. These Medical Devices (MD) combine data from novel sensors and existing modalities like scanners with elaborate software processing to assist caregivers in the same way Flight Management Systems help a pilot flying planes. For instance, Blue Ortho's MD allows performing Total Knee Arthroplasty (TKA, i.e. the placement of a femoral and tibial knee prosthesis) more precisely, potentially dividing the number of revisions by two [START_REF] Schnurr | Influence of computer navigation on TKA revision rates[END_REF]. Unfortunately, innovation and safety of such MD is hindered by the current software verification practices of the MD industry. Some formal methods have demonstrated their effectiveness for decades in the most safety-critical industries (defense, air, space, railways, nuclear power plants), but they failed to gain broader adoption.

MD provide an interesting opportunity for improving verification processes, not only during the manufacturers verification activities, but also in their clinical context for their whole lifetime. Indeed, these powerful systems can easily be equipped with elaborate tracing facilities, and exploiting their execution traces and sensor data can provide us with an unbiased and precise understanding of their behaviour in the field. Consequently, model-based verification of execution traces, also known as Runtime Verification, seems particularly fit to the task at hand. This strategy has, to our knowledge, never been applied to the MD industry, and a success in this highly visible domain could open the way to a broad range of formal methods and tools in numerous domains where only a small part of requirements are considered critical. For the MD industry, we believe runtime verification can advance the evaluation of MD safety and performance in the field. Although the generation of MD execution traces is a standard practice for a posteriori understanding when problems have been experienced, their systematic exploitation is new to the medical devices industry. A model-based, semi-automated analysis of these data provides a mean to verify that requirements are satisfied in any conditions, to validate the hypotheses made on the environment, and to understand the usage of the MD in a perspective of improvement. It can also bring disruptive progress in the postmarket surveillance of such systems which is usually limited to ineffective and informal followups of the device users. For instance, it will allow validating safety-related assertions on users behaviour [START_REF] Blandford | Who's looking? invisible problems with interactive medical devices[END_REF] or hardware longevity [START_REF] Sipos | Log-based predictive maintenance[END_REF]. Furthermore, some regulatory administrations like the FDA are interested in evidences of safety, i.e. factual data on the actual device.

The MODMED project gathers Blue Ortho, a MD manufacturer, MinMaxMedical, a software provider for MD and the LIG laboratory to propose a solution to modelbased verification of execution traces of MD in two folds: a structured trace library, designed with low adoption costs in mind, and an intuitive property specification language along with its implementation of verification tools.

Trace Library

Software execution is almost always traced to facilitate its development, however such traces are difficult to use for formal verification purposes because they are often incomplete and they lack abstraction.

Existing trace libraries such as ETW 1 , journald 2 , or log4j 3 are focused on performance or flexibility at runtime but very limited in structuring the data of logged messages. When not limited to ambiguous text or primitive data types, those libraries require more work than most developers can spend to format structured data and describe it in external data schemas. The most common trace attributes (severity, timestamps, etc.) are usually standardized as with syslog 4 , but there should also be an easy way to standardize the most prominent domain-specific data types and attributes (patient id, RMS, cut plane, medical terms, etc.) wherever they can be in the event data.

Few research is dedicated to improving trace generation tools and methodologies. Aspect-oriented programming tools like AspectJ [START_REF] Kiczales | An overview of AspectJ[END_REF] are fitted to automatically insert tracing code, but a survey of MinMaxMedical customers showed that medical software project managers are reluctant to use tools modifying the generated executable. They would prefer approaches like LogEnhancer [START_REF] Yuan | Improving software diagnosability via log enhancement[END_REF] that focus on improving traces quality before quantity. From their standpoint, the ideal solution would be to improve traces without even having to change trace points.

MODMED embraces this goal and proposes a C++ structured trace library that can greatly improve traces structure with existing trace points, and then, easily trace complex event data types without complicating trace points. This library 5 extends Qt logging 6 to:

• automatically collect and log event metadata such as category, function and abstract message format without requiring extra developer work or log space.

• structure event data using the source language trace expression and type system.

• unambiguously format event data in a TSV+JSON format fitted to human exploration and automated analysis (without having to maintain an external data schema).

• allow on-purpose classification of events based on their severity, category and abstract format message.

Expression and Evaluation of Properties

Trace verification can be automated when requirements are expressed as formal properties. Our main objective is to design a property specification language in such a way that it should be readable by quality engineers and easy to manipulate for software engineers when writing specifications. For this purpose, we analysed a large set of requirements for the aforementioned TKA to determine events of interests and how requirements rule the ordering of events. The result of our analysis is a language designed to concisely and intuitively capture the requirements for interventions and the desired behaviors of MD 7 .

The language is an adaptation of the specifications patterns proposed by Dwyer et. al. [START_REF] Dwyer | Patterns in property specifications for finite-state verification[END_REF]. Temporal properties can be expressed through a intuitive pairing of scope constraints and patterns. Scopes delimit the intervals of a trace where a property should hold, and patterns capture a requirement. Scopes may be freely nested and offer much more flexibility than the original specification patterns. Additionally, properties can be time-constrained through timed scopes, which not only delimits trace intervals according to events, but also according to elapsing time. Moreover, the language was built upon the idea that event parameters should have a central role. As a consequence, each time an event is referred to, its parameters may be extracted and constrained according to events seen previously. Let us consider an example of requirement for TKA and its formalisation in the proposed language:

Before starting an acquisition phase, all the necessary sensor types must have been detected. Although we believe that such properties are expressed in user-friendly manner, we are actively working on tools that will help express and understand them.

Conclusion

Prototype versions of the trace library and the language interpreter have been released under a free license 5 . The bridge between the two components is not automated yet, but they can be used as standalone tools. The verification of a few properties over a corpus of traces recorded during real surgeries already produced interesting results. For instance, it appears that one out of two surgeries includes at least one accidentally repeated click on the device's touchscreen, which impairs its usability. The frequentness of this defect, attributed to a hardware flaw in the device's stand, would not have been discovered without systematic and automatic verification of the execution traces.

  5. https://forge.imag.fr/projects/modmed/ 6. https://doc.qt.io/qt-5/qmessagelogger.html 7. https://sites.google.com/a/minmaxmedical.com/modmed/deliverables If several lists of required sensor types are present, only the last one before starting must be considered. The beginning of the acquisition phase, a new sensor detection and the list of required sensors are reflected through the events StartAcquisitions, SensorDetected and SearchSensors, respectively. The type of a sensor is given by the type parameter of the second event and the list of required trackers by the types parameter of the last one. The requirement can be formalized as follows:

	before each StartAcquisitions,
	given last SearchSensors ss,
	forall type in ss.types,
	occurrence_of SensorDetected sd
	where sd.type == type

Acknowledgments

This work is funded by the ANR project MODMED (ANR-15-CE25-0010).