Christina Konstantina

George Alexander Topka

François Chliavoras

Hugues Senocq

Diane Vergnes

Daniel Samelor

Constantin Sadowski

Brigitte Vahlas

Caussat

Topka

George Alexander Chliavoras

François Senocq

Hugues Vergnes

Diane Samelor

Daniel Sadowski

Constantin Vahlas

Brigitte Caussat
email: brigitte.caussat@ensiacet.fr

Large

Keywords: Chemical vapor deposition, Silicon dioxide, Ozone, TEOS, Computational Fluid Dynamics, Kinetic model

Introduction

The chemical vapor deposition (CVD) of SiO2 has been widely studied since thirty years [START_REF] Kotani | Low-temperature APCVD oxide using TEOS-ozone chemistry for multilevel interconnections[END_REF], mainly in semiconductor manufacturing as gate dielectrics [START_REF] Inoue | Low Temperature Growth of SiO2 Thin Film by Double-Excitation Photo-CVD[END_REF][START_REF] Ding | Low Dielectric Constant SiO2:C,F Films Prepared from Si(OC2H5)4/C4F8/Ar by Plasma-Enhanced CVD[END_REF], and more recently in membrane-based gas separation and energy harvesting [START_REF] Choi | Obtention of ZrO2-SiO2 hydrogen permselective membrane by chemical vapor deposition method[END_REF][START_REF] Liang | Experimental investigation on spectral splitting of photovoltaic/thermal hybrid system with two-axis sun tracking based on SiO2/TiO2 interference thin film[END_REF][START_REF] Azizi | Comparison of permeability performance of PEBAX-1074/TiO2 , PEBAX-1074/SiO2 and PEBAX-1074/Al2O3 nanocomposite membranes for CO2/CH4 separation[END_REF]. Other applications include the use of SiO2 in optoelectronics, or as a component in antibacterial food packaging films [START_REF] Varghese | Novel antibacterial silver-silica surface coatings prepared by chemical vapour deposition for infection control[END_REF] and food safety monitoring sensors [START_REF] Liu | Induced SERS activity in Ag@SiO2/Ag coreshell nanosphere arrays with tunable interior insulator: SERS activity in Ag@SiO2/Ag core-shell nanosphere[END_REF]. Moreover, core@SiO2-shell materials for solid phase extraction, catalyst support and SiO2-nanocomposites used for wastewater treatment are receiving increasing attention [START_REF] Zhang | Preparation of core-shell magnetic ion-imprinted polymer for selective extraction of Pb(II) from environmental samples[END_REF][START_REF] Klaigaew | Effect of preparation methods on activation of cobalt catalyst supported on silica fiber for Fischer-Tropsch synthesis[END_REF][START_REF] Zhang | β-SiAlON ceramic membranes modified with SiO2 nanoparticles with high rejection rate in oil-water emulsion separation[END_REF][START_REF] Li | Amino group functionalized SiO2@graphene oxide for efficient removal of Cu(II) from aqueous solutions[END_REF]. SiO2 films are optically transparent, smooth and dense, acting as good insulation layers [START_REF] Liu | Influence of SiO2 insulation layers thickness distribution on magnetic behaviors of Fe-Si@SiO2 soft magnetic composites[END_REF][START_REF] Fonrodona | Nanocrystalline top-gate thin film transistors deposited at low temperature by hot-wire CVD on glass[END_REF] and corrosion barriers [START_REF] Elam | Atmospheric pressure roll-to-roll plasma enhanced CVD of high quality silica-like bilayer encapsulation films[END_REF]. The tetraethyl orthosilicate (Si(OC2H5)4, TEOS) / ozone (O3) thermal CVD processing route presents marked advantages compared to the other existing pathways using TEOS such as TEOS pyrolysis, particularly in the ability to treat thermally sensitive substrates at temperatures lower than 300°C [START_REF] Juárez | Optical and electrical characterization of SiO2 films obtained by atmospheric pressure chemical vapor deposition[END_REF][START_REF] Okuyama | Gas-phase nucleation in the tetraethylorthosilicate (TEOS)/O3 APCVD process[END_REF].

Films deposited by thermal CVD present excellent and tunable step coverage [START_REF] Fujino | Step-Coverage Simulation for Tetraethoxysilane and Ozone Atmospheric Pressure Chemical Vapor Deposition[END_REF][START_REF] Yuan | Optimization of SiO2 film conformality in TEOS/O3 APCVD[END_REF][START_REF] Levin | The step coverage of undoped and phosphorus-doped SiO2 glass films[END_REF], while hydrogen content and porosity are lower than in films processed by PECVD [START_REF] Nagel | Porous SiO2 films prepared by remote plasma enhanced chemical vapour deposition -a novel antireflection coating technology for photovoltaic modules[END_REF], improving as such the material's barrier and insulation properties. Furthermore, in comparison with the thermal TEOS/O2 pathway, thermal TEOS/O3 leads to higher deposition rates at lower operating temperatures, thanks to the reactivity of O3 and the oxygen radicals generated through its thermal decomposition [START_REF] Okuyama | Gas-phase nucleation in the tetraethylorthosilicate (TEOS)/O3 APCVD process[END_REF]. The TEOS/O3 route is therefore promising to treat complex-inshape thermally sensitive substrates for innovative applications.

Succeeding in conformally coating complex substrates at low temperatures depends strongly on O3 concentration [START_REF] Shareef | Subatmospheric chemical vapor deposition ozone/TEOS process for SiO2 trench filling[END_REF][START_REF] Kawahara | Reaction mechanism of chemical vapor deposition using tetraethylorthosilicate and ozone at atmospheric pressure[END_REF]. The radicals produced by this chemistry are unstable, while the gas phase reactions are strongly dependent on temperature and radical species concentration [START_REF] Shareef | Role of gas phase reactions in subatmospheric chemical-vapor deposition ozone/TEOS processes for oxide deposition[END_REF]. Overcoming these drawbacks of the TEOS/O3 route requires good control over the interplay between the transport phenomena, i.e. the local gas velocity, temperature and reactive species concentrations, as well as the homogeneous and heterogeneous chemical reactions leading to the formation of the film. As this interplay is difficult to establish experimentally, a valuable approach is to associate a targeted experimental investigation of the growth rates with reactor modelling. In particular, Computational Fluid Dynamics (CFD) codes can be used, in which appropriate chemical reactions and kinetic laws are implemented. In this perspective, the identification of the involved reactions and kinetics is often the bottleneck, since a wealth of unknown homogeneous and heterogeneous reactions occur simultaneously in CVD processes. A classical way to overcome this difficulty is to consider only apparent chemical reactions, meaning that each reaction represents, in a simplified way, many real chemical reactions that are otherwise difficult to determine in detail. The advantage of such an approach is the relative rapidity of development and the accuracy of the kinetic model, while the main drawback is that it is only valid for the operating range, in particular that of the deposition temperature for which it has been developed.

For the TEOS/O3 CVD route, most of the relevant studies have been conducted in coldwall low or atmospheric pressure CVD reactors using showerhead injectors [START_REF] Schoof | Simulation and validation of SiO2 LPCVD from TEOS in a vertical 300 mm multi-wafer reactor[END_REF][START_REF] Adachi | Numerical Simulations of Films Formed by Cluster/Particle Co-Deposition in Atmospheric-Pressure Chemical Vapor Deposition Process Using Organic Silicon Vapors and Ozone Gas[END_REF][START_REF] Fujimoto | Effect of cluster/particle deposition on atmospheric pressure chemical vapor deposition of SiO2 from four gaseous organic Sicontaining precursors and ozone[END_REF], such as the industrial Watkins Johnson reactors treating 200 mm silicon wafers, which were the main reactors of interest for the microelectronic industry in the nineties. Kim and Gill were among the first authors to study the TEOS/O3 chemistry in a cold-wall reactor operating between 4 and 12 kPa, in a 280-405°C temperature range [START_REF] Kim | Modeling of CVD of Silicon Dioxide Using TEOS and Ozone in a Single-Wafer Reactor[END_REF]. They developed a global apparent mechanism in which TEOS decomposes in the gas phase in the fictitious molecule "INT" and in ethanal, through an excited O3 * centered mechanism. This INT molecule represents the various unknown chemical species produced by the homogeneous decomposition of TEOS. It can be formed through the substitution of TEOS' ethoxy ligands by hydroxyl ligands, following the formula Si(OH)n-1(OC2H5)5-n , where n = 1-4 [START_REF] Niéto | Analyse et modélisation d'un réacteur de CVD à fonctionnement continu : Application au dépôt de SiO2 pur et dopé[END_REF]. Additionally, INT can encompass oligomers, such as siloxane ones that are formed from the hydro-substitution of two intermediates with the previously mentioned formula [START_REF] Satake | Detection of intermediates in thermal chemical vapor deposition process using tetraethoxysilane[END_REF][START_REF] Sorita | The formation mechanism and step coverage quality of Tetraethylorthosilicate-SiO2 films studied by the micro/macro cavity method[END_REF]. A cascade of reactions resulting in the stepwise cleavage of the ethoxy ligands is also described [START_REF] Romet | Modeling of silicon dioxide chemical vapor deposition from tetraethoxysilane and ozone[END_REF]. Once formed, the INT species can either contribute to SiO2 deposition or decompose into parasitic gaseous by-products. Working between 400 and 425°C and for a pressure range of 10 3 -10 5 Pa, Dobkin et al. [START_REF] Dobkin | Mechanisms of deposition of SiO2 from TEOS and related organosilicon compounds and ozone[END_REF] modified this model by assuming that deposition occurs due to the presence of oxygen radicals, and that both TEOS and the INT species contribute to film formation. Zhou et al. [START_REF] Zhou | Numerical study of TEOS/O3 CVD mechanism in an industrial reactor[END_REF] developed a new model based on the two previous ones for a Watkins Johnson reactor operating at atmospheric pressure between 500 and 550°C. The main differences concern the assumption that TEOS must react with O3 to form, in one hand, the INT species and gaseous unreactive molecules noted R, and on the other hand the SiO2 film. The unreactive molecules R mentioned in the literature are CH3CHO, EtOH, MeOH [START_REF] Kawahara | Reaction mechanism of chemical vapor deposition using tetraethylorthosilicate and ozone at atmospheric pressure[END_REF][START_REF] Flores | Boundary layer chemistry probed by in situ infrared spectroscopy during SiO2 deposition at atmospheric pressure from tetraethylorthosilicate and ozone[END_REF][START_REF] Arnó | Fourier transform infrared characterization of downstream gas-phase species generated by tetraethylorthosilicate/ozone atmospheric pressure reactions[END_REF], while other authors assume the direct decomposition of such molecules into CO, CO2 and H2O [START_REF] Dobkin | Mechanisms of deposition of SiO2 from TEOS and related organosilicon compounds and ozone[END_REF][START_REF] Zhou | Numerical study of TEOS/O3 CVD mechanism in an industrial reactor[END_REF]. Nieto et al. [START_REF] Nieto | Modelling of an industrial moving belt chemical vapour deposition reactor forming SiO2 films[END_REF] applied this model to the same type of reactor and adapted some kinetic constants to its specific chamber configuration. In particular, they decreased the kinetics of the deposition reaction from TEOS, and the conversion of INT species into parasitic by-products that do not contribute to the deposition. It is worth noting that for these studies, no thermal activation appears for the kinetics of the deposition from TEOS. Romet et al. [START_REF] Romet | Modeling of silicon dioxide chemical vapor deposition from tetraethoxysilane and ozone[END_REF] improved the model of Dobkin et al. by considering up to 37 gas phase reactions, including the direct contribution of O3 to intermediate species formation, gas phase oligomerization and formation of diols. Similarly to previous studies, they did not consider any temperature dependence of the deposition reactions, and additionally removed any direct contribution of TEOS to the formation of the SiO2 film. However, their model failed to simulate the process above 400°C at low TEOS flow rates.

In Watkins Johnson type reactors, the distance between the gas injector and the wafer is lower than 20 mm and the residence time of the gaseous reactive phase is short (< 0.5 s), due to total flow rates exceeding 20000 standard cubic centimeters per minute (sccm). Moreover, the wall and substrate temperature ranges are spanning mainly between 400 and 550°C. It is thus inappropriate to use the developed chemical reactions and associated kinetic laws for significantly different configurations, in particular in reactors involving higher residence times and/or lower deposition temperatures, as we will present in the following sections.

The present work aims to develop a kinetic model for the atmospheric pressure CVD of SiO2 from TEOS, O2 and O3 mixtures, with a larger range of validity in terms of substrate temperature and residence time, more suited to CVD reactor geometries and deposition parameters able to treat large and/or complex-in-shape substrates. We apply a methodology combining experimental studies and process analysis by CFD simulation, in order to bring new information on the chemical and kinetic phenomena governing the deposition. In that which follows, the experimental and the computational details will be introduced first. Then, the experimental thermal profiles and associated deposition rates will be presented and discussed. These results will first be simulated using existing literature models. The aim is to reveal discrepancies, which will serve to develop a novel kinetic scheme. Finally, the elaborated kinetic model will be validated by comparison with the experimental information, and the simulation results will be discussed prior to providing concluding remarks.

Materials and methods

Experimental process

CVD of SiO2 films was performed in a horizontal, tubular, hot-wall reactor, schematically illustrated in Fig. 1a. The reactor consists of a fused silica tube with an inner diameter of 46 mm and a length of 700 mm, heated by a semitransparent furnace (Trans Temp, Thermcraft Inc.). The wall's temperature profile was measured every 20 mm along the tube under a N2 gas flow of 4036 sccm using a type-K thermocouple with an accuracy of ±2°C. Deposition experiments were performed at four different set point (SP) temperatures, namely 400, 450, 500 and 550°C, defined as the temperature measured at 360 mm from the inlet, after which a 140 mm long isothermal region exists. The CVD reactor inlet is connected to three gas lines. A N2 (99.9999%, Messer) dilution line with variable flow allows to fix the total gas flow at 4036 sccm for all runs and was heated at approximately 100°C. A second line bubbled 89, 223 or 444 sccm of N2 gas through a bubbler system containing TEOS (99.999%, Sigma-Aldrich), and transported the precursor's vapors towards the reactor. The bubbler was heated at 65.3±0.1°C, while the corresponding N2 line was kept at approximately 85°C to avoid recondensation of the precursor vapors. An O2 (99.9995%, Messer) line passed 1960 sccm through an ozonator (Lab2b Laboratory Ozone Generator, Triogen LTD). The ozonator was operating at its maximum ozone production, providing an O3 concentration of 60 mg/l at standard temperature and pressure, and thus supplying the reactor with an O2/O3 mixture. All three gas-lines were kept separated, and their contents were mixed just at the inlet of the reactor. The system was always in excess of O2. The reactor's base and operating pressures were monitored with a Pirani (MKS MicroPirani Transducer Series 925C) and a Baratron (MKS Baratron Type 627) gauges respectively, positioned downstream of the deposition zone, before a liquid nitrogen trap and a dry pump (Edwards, soXu 20iC). The operating pressure was fixed at 97.3 kPa (730 Torr) for all runs. 280 μm thick monocrystalline silicon (100) wafers (Neyco) were cut in rectangles of 32×24 mm 2 and used as substrates. Before being inserted into the reactor, the substrates were degreased in a succession of three ultrasound bath steps, including: 1) an ultrasound bath using distilled water for 5 min, then rinsed with acetone, 2) an ultrasound bath using acetone (>99%, VWR Chemicals) for 5 min, then rinsed with ethanol, 3) an ultrasound bath using ethanol (99.3%, VWR Chemicals) for 5 min, and finally dried under Ar (99.9999%, Messer) flow.

Within the reactor, the substrates were supported vertically by home-made, tubular, stainless steel substrate holders, with a 7 mm deep insertion slot, positioned on the row. Their total length was 450 mm, supporting a maximum of 18 coupons for each run. The first substrate was placed at a distance of 57 mm from the reactor inlet. Similar to Ponton et al.'s work [START_REF] Ponton | Development of a kinetic model for the moderate temperature chemical vapor deposition of SiO 2 films from tetraethyl orthosilicate and oxygen[END_REF], the goal was to process SiO2 films deliberately in a 360 mm long non-isothermal zone, existing before the ca.140 mm long isothermal zone, in order to obtain complementary mechanistic and kinetic information. Moreover, placing the first few substrates close to the reactor inlet allows addressing and quantifying the direct contribution of TEOS to the deposition, since the TEOS concentration is known only at the reactor inlet. Six runs were performed with a deposition duration ranging from 30 to 90 min, as summarized in Table 1. Through experiments E1 to E3, the influence of the TEOS flow rate was evaluated for a SP temperature of 400°C. Through runs E3 to E6, the effect of the SP temperature was investigated, by fixing all flow rates.

Experimental Validation

The average SiO2 deposition rates were calculated from the deposited mass, by measuring each sample's weight difference before and after deposition using a Sartorius Genius Series ME215P analytical balance with 0.015 mg precision. For simplicity, it was assumed that deposition on each substrate surface is uniform. To convert the deposited mass to thickness and deposition rate, the films' density was taken equal to 2200 kg.m -3 [START_REF] Nieto | Modelling of an industrial moving belt chemical vapour deposition reactor forming SiO2 films[END_REF]. More accurate local thickness profiles were measured by reflective ellipsometry, using a Semilab SE-2000 ellipsometer operating in the 250-1000 nm wavelength range. The incidence angle was kept constant at 70°. Each sample was probed in the center, 2 mm below the top edge, and the spectroscopic ellipsometry data were fitted in the 250-1000 nm wavelength range using the software Semilab SEA. The Sellmeier model was utilized to extract the measured local thicknesses values. Fourier Transform Infra-Red (FTIR) spectroscopy measurements, detailed in Fig. S1 of Supplementary Materials, revealed that for all the conditions tested, the films are composed of pure SiO2. The total reactive volume of the CVD reactor was represented by a three-dimensional (3D) geometrical domain (Fig. 1b) of 1,927,954 hexahedral and tetrahedral cells. Due to symmetry along the YZ plane (top right inset scheme of Fig. 1) and in order to reduce the computational time, only half of the reactor was considered.

Process model

A
The following boundary conditions were applied:

-A flat profile was considered for the gas velocity at the gas inlet. The total mass flow rate, and all species' mass fractions were fixed to the experimental values. The inlet gas temperature was set equal to the temperature of the wall close to the inlet.

-A symmetry boundary condition was applied at the YZ symmetry plane.

-A classical no-slip condition was imposed for the gas velocity on the reactor walls, the substrate holder, and the substrates themselves. The temperature of the surfaces was set equal to the temperature profile measured experimentally on the reactor walls before deposition. Due to the steep decrease of the temperature at the exhaust, a constant temperature was applied for the last 100 mm of the reactor; i.e. between 600 and 700 mm, to avoid backflow in the computations and to speed up convergence. Modifying the temperature profile in this region does not impact the results, because samples were placed only between 57 and 500 mm from the reactor inlet. The mass flux density of each species was assumed to be equal to the corresponding heterogeneous reaction rate.

-The reactor exhaust was defined with a pressure-outlet boundary condition. The total pressure was fixed at the operating pressure. A zero diffusion flux was fixed for all gas phase variables.

The chemical reaction rates were implemented into FLUENT ® in the form of an Arrhenius type expression:

Rr=𝑘 0,𝑟 𝑒𝑥𝑝 (- 𝐸 𝑎 𝑟 𝑅 𝑇 ) ∏ 𝑝 𝑖 𝑖=𝑛 𝑖=1 (1)
For the surface reaction rates, Rr (kmol.m -2 .s -1 ) stands for the rate of the reaction r, k0

(kmol.m -2 .s -1 .Pa -n ) is the pre-exponential factor of the reaction r, Ea r (kJ.mol -1 ) is the activation energy of the reaction r, T (K) is the temperature of the surface where the reaction r occurs, n is the number of involved species, pi (Pa) is the concentration of the species i. For the homogeneous chemical reactions, Rr is expressed in kmol.m -3 .s -1 and k0,r in kmol.m -3 .s -1 .Pa -n , whereas T (K) is the temperature of the gas phase.

Results

Experimental thermal profiles and deposition rates

Fig. 2 shows photographs of the SiO2 films deposited across the first 12 Si samples placed in the CVD reactor, for the four SP temperatures tested, corresponding to runs E3 to E6.

The deposition rate on substrates 13 to 18 was similar to the one on substrate 12. No powder was present on the films. Radial and axial iridescences appear on the surface of the substrates, indicating thickness variations of the films. These variations are directly linked to the chemical mechanisms responsible for the film deposition, as it will be detailed in the forthcoming sections, and are probably enhanced by the thermal gradients existing in the first part of the reactor. The corresponding simulated thickness profiles are also presented in Fig. 2 and will be discussed in section 3.4. The thermal profiles and the experimental deposition rates are plotted in Fig. 3 versus the distance from the reactor inlet, for runs E3 to E6. For all four runs, the temperature in the first 70 mm of the reactor remains low, i.e. between 115 and 156°C. After 70 mm, the temperature increases rapidly to reach the target SP temperature at 360 mm. Then, an isothermal zone appears, spanning approximately from 360 to 500 mm. For run E3, the deposition rate first increases to reach 9 nm.min -1 until the temperature reaches approximately 250°C, then decreases in the zone where temperature continuously increases. For runs E4 to E6, the deposition rate reaches a maximum value of 11-13 nm.min -1 at approximately 180°C, after which it decreases with the temperature, and remains constant within the isothermal region. This behavior could correspond to systems for which at least two different species are responsible for deposition, as it will be discussed in the forthcoming sections.

Simulation of the deposition rate using existing model

Nieto et al.'s model [START_REF] Nieto | Modelling of an industrial moving belt chemical vapour deposition reactor forming SiO2 films[END_REF] was used as the starting point for the present work, because, to the best of our knowledge, it is the study closest to the conditions used in our experiments. Their kinetic system (Table S1 in the Supplementary materials) was implemented into FLUENT ® and run E3 was simulated (Fig. 4). A high deposition rate (33 nm.min

Development of a new kinetic scheme

According to the previous preliminary experimental observations and computational results, we modified Nieto et al.'s model and kinetic laws based on the following assumptions:

(i) the deposition in the vicinity of the inlet region is attributed to the reaction between TEOS and O3 (R6), exclusively; (ii) the deposition closer to the exhaust, where O3 is no more present, is attributed to INT.

First, concerning the kinetics of (R6), the experimental results of runs E1 to E3 were considered in order to calculate a new partial order exponent a for TEOS. From each run, performed at the same temperature, the first sample from the inlet was selected and its weight gain after deposition was converted into deposition rate. Runs E1 to E3 are coupled in pairs of two, and by using the (R6) kinetic law of Table S1 (Supplementary Materials), the ratio of their reaction rates is expressed as:

𝑅6 𝐸 𝑖 𝑅6 𝐸 𝑗 = 𝑘 6 [𝑇𝐸𝑂𝑆] 𝑖 𝑎 [𝑂 3 ] 𝑖 0.25 𝑘 6 [𝑇𝐸𝑂𝑆] 𝑗 𝑎 [𝑂 3 ] 𝑗 0.25 (2)
where R6 is the reaction rate of reaction (R6), Ei and Ej are the experimental runs E1, E2, E3, with i = 1-2, j = 2-3, and i ≠ j. k6 is the kinetic constant of reaction (R6), [TEOS]i and

[TEOS]j are the concentrations of TEOS at the inlet for runs Ei and Ej respectively, [O3]i and

[O3]j are the concentrations of O3 at the inlet for runs Ei and Ej respectively, and a is the partial order exponent of TEOS.

With the kinetic constant k6 remaining unchanged (constant temperature for the considered runs), and assuming that the concentration of O3 at the inlet is approximately equal across all runs, equation ( 2) is simplified to equation (3).

𝑅6 𝐸 𝑖 𝑅6 𝐸 𝑗 = ( [𝑇𝐸𝑂𝑆] 𝑖 [𝑇𝐸𝑂𝑆] 𝑗 ) 𝑎 (3)
Equation ( 3) is then solved for a. The logarithm of

𝑅6 𝐸 𝑖 𝑅6 𝐸 𝑗 is plotted versus the logarithm of [𝑇𝐸𝑂𝑆] 𝑖
[𝑇𝐸𝑂𝑆] 𝑗 (Fig. 5a). From the slope of Fig. 5a, we deduce a value for a equal to 1.9, with reaction (R6) now becoming:

𝑅6 = 𝑘 6 [𝑇𝐸𝑂𝑆] 1.9 [𝑂 3 ] 0.25 (4)
This value is higher than the one used by Our experimental data show a temperature-dependent evolution of the deposition rate in the region where (R6) is active. This can be observed in Fig. 3 when comparing the deposition rates of the first sample for each experimental run, a location where the partial pressures of TEOS and O3 in the gas phase are assumed approximately similar across all runs, the temperature being the only variable. Therefore, the first and second samples closest to the inlet were selected from runs E3 to E6, and their respective deposition rate was divided by

[𝑇𝐸𝑂𝑆] 𝑠 1.9 [𝑂 3 ] 𝑠 0.25 (by solving equation ( 4) for k6) and plotted versus the inverse of their local average temperature, spanning from 115 to 230°C (Fig. 5b). It is worth noting that the two points relative to the lowest SP temperature experiment (run E3) are not aligned with those of the other runs. It is possible that the mechanisms involved at temperatures lower than 150°C are different than those prevailing between 150°C and 230°C. Thus, these two points were not considered. As a consequence, the final model will have a lower accuracy for deposition temperatures lower than 150°C. The temperature dependent kinetic rate for (R6) was obtained from the results of runs E4 to E6, for which an activation energy of 21.8 kJ.mol -1 was deduced.

It is difficult to comment on this value in relation to literature data, because in this kinetic model, each apparent chemical reaction represents a set of multiple real chemical reactions.

Moreover, it is the first time that an activation energy has been calculated for this reaction. Unlike Romet et al. [START_REF] Romet | Modeling of silicon dioxide chemical vapor deposition from tetraethoxysilane and ozone[END_REF] who proposed a deposition mechanism involving the reaction of silanol or diol intermediates with O3 and O • , we used a reaction between INT with O2 instead, as it represents our experimental results more accurately. From the above, an activation energy of 22.4 kJ.mol -1 was deduced, and reaction (R7) was reformulated, as detailed in Table 2. Since no decline in the deposition rates was observed close to the exhaust, reaction (R5) was removed and not applied in our model, similarly to the model of Kim and Gill [START_REF] Kim | Analytical model for chemical vapor deposition of SiO2 films using tetraethoxysilane and ozone[END_REF]. The final version of our model and the respective kinetic constants are given in Table 2.

It is recalled that all the chemical reactions considered in the model are apparent and stand for numerous elementary mechanisms. Phenomena such as reactive species adsorption, surface diffusion, chemical reaction and product desorption are accounted for in one single surface reaction.

Table 2. Proposed apparent chemical model, with the respective kinetic laws and constants.

𝐑𝐞𝐚𝐜𝐭𝐢𝐨𝐧 𝐒𝐭𝐨𝐢𝐜𝐡𝐢𝐨𝐦𝐞𝐭𝐫𝐲 𝐑𝐞𝐚𝐜𝐭𝐢𝐨𝐧𝐑𝐚𝐭𝐞 𝐚𝐧𝐝 𝐊𝐢𝐧𝐞𝐭𝐢𝐜 𝐜𝐨𝐧𝐬𝐭𝐚𝐧𝐭𝐬

Kinetic constants units

Volumetric

Reactions

(𝑘𝑚𝑜𝑙 𝑚 3 𝑠 ⁄ ) (R1) 𝑂 3 + 𝑀 → 𝑂 2 + 𝑂 • + 𝑀 𝑘 1 [𝑂 3 ][𝑀] 𝑘 1 = 2.5 × 10 11 × 𝑒𝑥𝑝 (- 11430 𝑇(𝐾) ) 𝑚 3 𝑘𝑚𝑜𝑙 -1 𝑠 -1 (R2) 𝑂 3 + 𝑂 • → 2𝑂 2 𝑘 2 [𝑂 3 ][𝑂 • ] 𝑘 2 = 10 10 × 𝑒𝑥𝑝 (- 2090 𝑇(𝐾) ) 𝑚 3 𝑘𝑚𝑜𝑙 -1 𝑠 -1 (R3) 2𝑂 • + 𝑀 → 𝑂 2 + 𝑀 𝑘 3 [𝑂 • ] 2 [𝑀] 𝑘 3 = 4 × 10 8 × 𝑒𝑥𝑝 (+ 720 𝑇(𝐾) ) 𝑚 6 𝑘𝑚𝑜𝑙 -2 𝑠 -1 (

Validation of the proposed kinetic model

Fig. 7 presents the calculated and experimental deposition rate profiles along the reactor length with the respective TEOS and INT contribution to the film for runs E3 to E6. It appears that the new kinetic model provides very good agreement with the experimental data, with the average relative error of all runs being approximately 17% for the samples located in the isothermal region. The poorest agreement is observed at the inlet zone of run E3 (Fig. 7a) corresponding to the lowest SP temperature, i.e. 400°C. This is not surprising since we did not consider the inlet zone results of this run to build our model, as explained in section 3.3. The overestimation of the experimental results in this region further underlines that at a different mechanism possibly takes place for temperatures lower than 150°C (Fig. 5b).

For runs E4 to E6, (Fig. 7b to 7d), the agreement in the inlet zone is clearly improved with increasing SP temperature. For run E5, i.e. SP 500 °C, experimental errors appear mainly in the isothermal zone due to a slight vertical angular displacement of specific substrates within the substrate holder, causing them to be out of vertical alignment. Lastly, for run E6, i.e. 550°C (Fig. 7d), a shoulder is observed on the experimental deposition rate data, spanning between 150 and 275 mm, corresponding to a zone where the temperature reaches 520°C. The model fails to conform to this deposition profile, underestimating the total deposition rate in this region. The existence of this shoulder cannot be attributed to experimental error, such as the vertical, out-of-alignment displacement of the samples, as it is rather uniform in shape and does not appear to be random. It is likely the result of another intermediate reaction contributing to deposition, activated starting from 520°C, meaning that the model is less accurate above this temperature.

Observing TEOS' and INT's respective contributions to SiO2 deposition, the simulation results confirm that TEOS is responsible for the formation of the film in the coldest zones close A remarkable agreement is observed between the calculated and experimental results of Fig. 8, confirming the ability of the model to predict local values of deposition rates. For the first two coupons, the simulations overestimate the local thickness due to the low temperature, as previously explained. Despite that, Fig. 8 is useful for ascertaining that the model can give

The W shape is due to the two-portion split of the gas previously mentioned and is also observed for the mass fraction profiles as detailed in the Supplementary materials. The bulk of the gas is logically colder than the gas close to the reactor walls and the substrates. Such profiles can be explained by the dominance of the convective phenomena and by the low thermal conductivity of the gas. Fig. 9 reveals that the O3 and TEOS mass fractions remain constant until the average gas temperature reaches a value close to 200°C at approximately 140 mm from the inlet. From this zone onwards, the TEOS mass fraction decreases due to the heterogeneous reaction (R6), and then a little further into the reactor, due to the homogeneous (R4) reaction (Fig. 10). O3 is depleted slightly more rapidly in the same zone, until only a very low mass fraction remains. This depletion is attributed mainly to reaction (R1), whose kinetics is much higher than all the other homogeneous reactions (Fig. 10), and also to reactions (R4) and (R8). Reaction (R4) gives rise to the INT species formation in a localized region between approximately 120 and 240 mm from the inlet, after which the reaction rate of (R4) diminishes to very low value, because no more O3 is present after this zone. For the same reason, reaction (R8) is also active only in this region, as a priori defined. After about 280 mm from the inlet, the TEOS and O3 mass fractions remain constant, with the O3 mass fraction being close to zero, and as such reactions (R4) and (R6) are unable to take place. On the contrary, INT is slowly depleted through reaction (R7), and presents a remarkable 200 mm long zone of uniform value, between 250 and 450 mm from the inlet. Interestingly, despite the mass fractions of INT and O2 remaining fairly constant after 450 mm, the reaction rate of (R7) decreases, underlining its temperature dependency. The atomic oxygen is formed through reaction (R1) just after the zone corresponding to the concentration decrease of O3 (between approximately 120 and 220 mm), and appears in low concentrations for most of the reactor's length. For O2, it is present in excess and thus its mass fraction is uniform inside the reactor, with a slight increase observed in the O3 depletion zone, mainly due to reaction (R1) and also to reactions (R2) and (R3).

Overall, the TEOS conversion rate for run E4 is approximately 79%, increasing for higher SP temperatures. Despite that, the global TEOS-to-SiO2 conversion rate on all solid surfaces remains low does not exceeding 5.2% for all runs. This means that the TEOS conversion involves mainly the gas phase reaction (R4), which converts TEOS into INT. The key information is that reaction (R7) should be favored to uniformly coat complexin-shape and/or large substrates. In the frame of the development of a new process, the model can be used to optimize both the CVD reactor design and the deposition conditions to enhance this heterogeneous reaction and decrease the other ones.

It is also worth noting the rapid consumption of O3 even at low temperatures. As a result, TEOS and INT remain unreacted, with reactions (R4), (R6) and (R8) ending prematurely. In those reactions, O3 depletion becomes the rate-limiting step, indicating that a higher O3 flow would be preferable all across the substrate zone. However, increasing the O3 concentration would lead to higher deposition rates across the whole reactor length in general, and especially at the inlet region where reaction (R6) rapidly takes place. As such, big discrepancies in deposition rates between the colder inlet zone and the isothermal region would still exist. An alternative studied was the increase of the total flow rate, as shown in Figure S7 of Supplementary Materials. By increasing the total flow rate up to three times, the deposition profile in the inlet zone is flattened and spread out over a larger region, thus demonstrating the capability of the kinetic model to optimize the deposition process. Nevertheless, the increase of the temperature as the gas progresses towards the isothermal region still leads to the thermal decomposition of O3 and its depletion in the later parts of the reactor. It becomes apparent that the combined contribution of reactions (R6) and (R7), along with reaction (R4), is required for more flat profiles. A proposed solution for higher TEOS-to-INT conversion rate, and TEOS and INT conversion into SiO2 film to result in higher deposition rates in the isothermal region and an overall flatter deposition profile would be a sequential introduction of ozone at multiple parts along the reactor. For better uniformity, it is proposed to insert O3 in already heated zones, in order to avoid the abrupt deposition rate due to (R6) in cold zones. This is a known problem, and such a solution has been proposed in other works [START_REF] Yuan | Optimization of SiO2 film conformality in TEOS/O3 APCVD[END_REF][START_REF] Nieto | Modelling of an industrial moving belt chemical vapour deposition reactor forming SiO2 films[END_REF]. Since a multiple injection setup of O3 cannot be easily implemented experimentally, the present model could be used to help finding the best geometry and operating conditions to increase intermediate conversion and overall yield.

Conclusions

The deposition of SiO2 from a TEOS/O2/O3 gas mixture performed by atmospheric pressure CVD is a promising route to provide new surface properties to thermally sensitive substrates. An apparent kinetic model is proposed, in order to represent CVD reactor geometries and deposition conditions able to treat such substrates. The range of validity of the new model is larger in substrate temperature (150-520°C) and in gas residence time (>0.5 s) than those of the literature.

We applied a methodology combining experimental studies and process analysis by CFD simulation. Experimental results obtained in non-isothermal zones of a hot-wall tubular reactor were simulated using the CFD code FLUENT ® .

Based on the model of Nieto et al. [START_REF] Nieto | Modelling of an industrial moving belt chemical vapour deposition reactor forming SiO2 films[END_REF], a thermal dependency of the direct formation of SiO2 from TEOS and O3 was considered for the first time. For the conditions tested, it was found that the intermediate species resulting from the gas phase decomposition of TEOS participate in two different apparent surface reactions, both producing SiO2 films. An agreement in shape and value of the calculated local profiles of deposition rates with experimental values validated the new kinetic model.

The analysis of the calculated local evolutions of gas velocity, temperature and species mass fractions provides valuable information about the process behavior. In particular, the detailed contribution of each species to the deposit was discussed and alternative routes of process optimization were proposed. It appears that the model can help in determining the best reactor design and operating conditions to uniformly coat complex and/or large substrates, even at atmospheric pressure. The model could also help in increasing the deposition rate and the TEOS contribution to the deposit by designing a multiple injection system of ozone. This new model paves the way to develop and optimize new CVD processes allowing to coat a large variety of substrates by SiO2 thin films at quite low temperatures for innovative applications.

Fig. 1 .

 1 Fig. 1. Schematic (a) and 3D (b) representation of the CVD reactor (not true-to-size and

  CFD model considering gas phase and surface chemical reactions was developed to represent the reactive transport phenomena and the deposition rate taking place in the CVD reactor. The local gas flow and temperature profiles in the reactor were calculated by solving the mass and momentum conservation and thermal energy balance equations in each elementary discretized volume. The local distribution of species mass fractions and SiO2 deposition rate were calculated after considering the kinetic data of the homogeneous and heterogeneous reactions into the species conservation equations. The implemented chemical reactions and kinetic laws were based on those developed by Zhou et al.[START_REF] Zhou | Numerical study of TEOS/O3 CVD mechanism in an industrial reactor[END_REF] and Nieto et al.[START_REF] Nieto | Modelling of an industrial moving belt chemical vapour deposition reactor forming SiO2 films[END_REF], as detailed below. The numerical simulations were run using ANSYS ® FLUENT ® 18.2, with a cellcentered finite volume scheme and second-order spatial accuracy. The assumptions considered to simplify the numerical complexity are the following:-Steady state regime, -Laminar gas flow (Reynolds number lower than 1,000), -Incompressible gas flow, due to the low value of the Mach number (<0.04) in the reactive zone, -Heat of reactions neglected, due to high N2 dilution, -Ideal gases.The physical properties of O2, O3, H2O, CO, CH3CHO and N2 were taken from the FLUENT ® database. They are calculated using the kinetic theory of gases and the Chapman-Enskog theory, as detailed in Bird et al.[START_REF] Bird | Transport Phenomena[END_REF], in which the Lennard Jones (LJ) parameters of these classical species are also given. The LJ parameters of the oxygen radical O • (σ = 2.75 Å and ε/k = 80 K), and of TEOS (σ = 7.03 Å and ε/k = 522.7 K) were taken from the CHEMKIN database[START_REF] Kee | A Fortran computer code package for the evaluation of gas-phase, multicomponent transport properties[END_REF]. Those of the intermediate (INT) were set equal to the parameters of TEOS. The molecular weight of INT was fixed at 180.27 g.mol -1 by matching it to triethoxysilanol, ((Si(OH)(OC2H5)3), one of the most probable intermediate species as reported by many authors[23,,32,35,41]. All physical properties of the by-product produced by the INT were set equal to the latter.

Fig. 2 .

 2 Fig. 2. Photographs of the SiO2 films along the substrate holder for runs E3 to E6, and

Fig. 3 .

 3 Fig. 3. Evolution of the experimental deposition rates (blue circles), the reactor's thermal

- 1 )

 1 was observed in the first 150 mm of the reactor, due to the intense surface reaction R6 (TEOS+6O3→SiO2+O2+byproducts) between TEOS and O3. No such deposition rate values were measured experimentally, as shown in Fig. 3a. This considerable discrepancy can be explained by the fact that Nieto et al. did not consider deposition on the reactor walls, and the average gas temperature of the single-wafer treated in Nieto et al.'s reactor was always about 400°C, a temperature at which the concentration of O3 is decreased due to its thermal decomposition through reactions R1 (O3+M→O2+O • +M), R2 (O3+O • →2O2) and R3 (2O • +M→O2+M). In contrast, the temperature of the first silicon substrate in run E3 does not exceed 150°C. Because k6 does not consider any temperature influence in Nieto et al.'s model, their kinetic mechanism is unable to represent our experimental results.

Fig. 4 .

 4 Fig. 4. Average calculated deposition rates of run E3 (SP 400°C), simulated using

  Zhou et al. (a = 1) or Nieto et al. (a = 0.4). The difference may be due to the steep thermal gradient observed in their reactor configurations, which activates the conversion of TEOS to INT through reaction (R4) more rapidly, downplaying the role of reaction (R6) and subsequently the influence of the TEOS concentration on the direct SiO2 deposition.

Fig. 5 .

 5 Fig. 5. (a) Determination of partial order exponent for TEOS in reaction (R6). (b) Linear

Fig. 6 Fig. 6 .

 66 Fig.6illustrates the experimental and the computed deposition rates along the length of

Fig. 7 .Fig. 8 .

 78 Fig. 7. Comparison of average experimental deposition rates (black squares) and average total

Fig. 9 .

 9 Fig. 9. Local velocity, temperature and mass fraction profiles of TEOS, INT, O3, O2 and O •

Fig. 10 .

 10 Fig. 10. Reaction rate profiles for all reactions plotted along the XZ plane for the gas phase

Table 1 .

 1 Deposition conditions of experimental runs performed in thermal CVD reactor.

	Run SP
	T
	(°C)

N2 flow rate through TEOS bubbler (sccm) TEOS flow rate (sccm) N2 dilution flow rate (sccm)

						Deposition
						time
						(min)
	E1	400	444	10	1632	30
	E2	400	223	5	1853	30
	E3	400	89	2	1987	50
	E4	450	89	2	1987	90
	E5	500	89	2	1987	30
	E6	550	89	2	1987	30

Acknowledgements

The present work was funded by ANR (Agence Nationale de la Recherche) under the contract HEALTHYGLASS ANR-17-CE08-0056. The authors gratefully acknowledge ANR for supporting the present research. The authors are also thankful to Irea Touche from LGC for her help on the numerical frame and troubleshooting, to Nikolaos Marios Katritsis for his aid in the initial implementation of the kinetic model, and to Olivier Debieu for his guidance on the ellipsometric measurements.

agreeable local thickness values, complementary to the average deposition rates calculated per coupon in Fig. 7. Lastly, the good agreement in both shape and value thus validates the new chemical and kinetic model as an appropriate simplified representation of the involved chemical mechanisms, for the range of operating conditions tested.

Discussion on process simulation results

The validated model provides valuable local hydrodynamic, thermal and mass transfer information for the studied chemical system. Taking run E4 as an example, Fig. 9 The velocity profile presents a double vortex-like shape in the inlet region and across the first substrates, because the gas is concomitantly heated by the solid surface of the wall and the substrates, and is split in two portions when arriving near the substrates (results detailed in Supplementary Materials). When plotting using the horizontal plane view of Fig. 1, the velocity profile appears parabolic, as classically observed for laminar gas flow. The hottest zones of the reactor logically give rise to the maximum velocity, close to 0.5 m.s -1 . The corresponding minimum residence time for the gas into the reactor is about 1.4 s. This is a high residence time for a hot wall CVD reactor, and could favor homogeneous chemical reactions, already enhanced by the operation at nearly atmospheric pressure.

The gas temperature remains below 150°C before the substrates' zone. It then increases and presents a W shape, followed by a parabolic shape in the zone occupied by the substrates.


This special issue contains a selection of papers originally presented at the 17th and 18th editions of the International Conference on Hybrid Systems: Computation and Control (HSCC'14 and HSCC'15).

HSCC has long been a leading, single-track conference on foundations, techniques, and tools for analysis, control, synthesis, implementation, and applications of dynamical systems that exhibit continuous and discrete (hybrid) dynamics. Applications include cyber-physical systems (CPS), mixed signal circuits, robotics, large-scale infrastructure networks, as well as natural systems such as biochemical and physiological models. In particular, HSCC solicits theoretical as well as applied research papers that present original work combining ideas from computer science and control systems.

HSCC'14 was held in Berlin, April 14-17, 2014, as part of the 7th CP-SWeek (Cyber Physical Systems Week) that in addition to HSCC comprised the Real-Time and Embedded Technology and Applications Symposium (RTAS), the International Conference on Information Processing in Sensor Networks (IPSN), the International Conference on Cyber-Physical Systems (ICCPS) and the Conference on High Confidence Networked Systems (HiCoNS). About 100 participants registered for the HSCC conference. The program of the conference, chaired by Martin Fränzle and John Lygeros, consisted of the presentation of 29 regular papers, 2 tool papers, 1 case study paper, 1 keynote lecture by Gregory Batt of the Institut National de Recherche en Informatique et Automatique (INRIA) on the topic of "Cells driven by computers: long-term model predictive control of gene expression in yeast" and a panel discussion on "The Role of Robustness in Hybrid Systems". Additionally, a total of 13 posters and demos were accepted by HSCC'14 for presentation at the joint CPSWeek poster and demo session.

HSCC'15 was held in Seattle, April 13-16, 2015, alongside RTAS, IPSN and ICCPS, as part of the 8th CPSWeek. About 120 participants registered

for the HSCC conference. The program of the conference, chaired by Antoine Girard and Sriram Sankaranarayanan, consisted of the presentation of 26 regular papers, 2 tool papers and 2 case study papers, 1 keynote lecture by Jyotirmoy V. Deshmukh from Toyota on the question"Will future cars have formally verified powertrain control software?". Additionally, a total of 11 posters and 6 demos were accepted by HSCC'15 for presentation at the joint CPSWeek poster and demo session.

The selection of papers for this special issue was made based on the reviewer evaluations of the conference papers, on the guest editors assessment and with the coverage objective of reflecting the diversity of the conference. The authors were asked to extend their papers to include new material with respect to the conference paper. These extended papers went through the standard review process of Nonlinear Analysis: Hybrid Systems, handled by the guest editors. At the end of the review process, ten papers were accepted for publication in the special issue. They cover the whole range of topics that touch upon hybrid systems, bridging control theory on the one end (for example, stability conditions for switched and networked control systems) to theoretical computer science on the other (for example, abstraction, verification and synthesis of hybrid systems). The issue also includes a paper on an exciting application of hybrid systems theory to bipedal locomotion. A short overview of the individual papers of the special issue is included below.

Papers in special issue

2.1. "Deciding the boundedness and dead-beat stability of constrained switching systems" by M. Philippe, G. Millerioux and R.M. Jungers In this paper, the authors investigate two problems related to stability analysis of discrete-time linear switched systems where switching signals are constrained by some finite transition graph. Sufficient conditions are established for boundedness of trajectories of constrained switched systems, when the constrained joint spectral radius of the system matrices is equal to one. An algorithm to check these conditions is provided. Then, dead-beat stability (i.e. finite time convergence to zero) is investigated and a polynomial time algorithm is provided for verifying this property.

"

Computing compositional proofs of Input-to-Output Stability using SOS optimization and δ-decidability" by A. Murthy, Md. A. Islam, S.A. Smolka and R. Grosu In this paper, the authors compute bisimulation functions to characterize input-to-output stability of systems through the combined use of Sum-Of-Squares programming and delta-decision procedures. The resulting bisimulation functions are used in conjunction with a small-gain theorem to replace components in a feedback composed system using their bounded error abstraction for the purposes of verification. The authors apply their framework on a bounded error cardiac myocyte model to replace a four variable model safely by a new model with a single variable, while proving bounded error guarantees. Extensions to non-polynomial vector fields are also presented.

"Finite abstraction of mixed monotone systems with discrete and con-

tinuous inputs" by S. Coogan and M. Arcak In this paper, the authors propose an efficient abstraction technique for a class of discrete-time nonlinear systems. The considered class is that of mixed monotone systems, an extension of the class of usual monotone systems. For such systems, the reachable set from a rectangular set can be over-approximated by evaluating a decomposition function at only two points. This property makes it possible to compute discrete abstractions of mixed monotone systems in a simple and efficient way. Both the case of discrete and continuous inputs are considered. Case studies in insect population dynamics or vehicular traffic networks are reported.

2.4. "Temporal logic control for stochastic linear systems using abstraction refinement of probabilistic games" by M. Svoreňová, J. Křetínský, M. Chmelík, K. Chatterjee, I. Černá and C. Belta In this paper, the authors consider the control of discrete-time, stochastic linear dynamical systems for temporal logic objectives written in the generalized reactivity (GR(1)) class of formulas. They provide an iterative approach that partially computes the set of winning initial states for the specification after any finite number of iterations. Further, upon termination, their procedure is guaranteed to yield the full set. The complexity of the algorithm is characterized. For the special case of reachability properties, a more efficient algorithm is provided. The algorithms are demonstrated on an illustrative case study.

"Bounded invariant verification for time-delayed nonlinear networked

dynamical systems" by Z. Huang, C. Fan and S. Mitra This article deals with the problem of verifying bounded invariant for dynamical systems featuring delayed interconnections, as found for example in networked control systems. It employs a notion of input-to-state discrepancy to generate a low dimensional system whose trajectories are bounding the distance between the trajectories of a complex network that may arise under changed initial states. This facilitates a sound and relatively complete algorithm for bounded invariant verification of networked nonlinear dynamical systems interacting through possibly delayed signals.

2.6. "On stability of discrete-time switched systems" by A. Kundu and D. Chatterjee This paper deals with stability conditions, primarily for switched nonlinear, discrete time systems. The discrete time switched systems considered are such that any point in time their state evolves according to dynamics selected among a finite family of transition functions, with the switching among the different possible dynamics controlled through a so-called switching signal. The authors first establish conditions on the switching signals that ensure global asymptotic stability to the zero state. In the spirit of dwell time results, the conditions are based on growth relations among multiple Lyapunov functions; later in the paper these conditions are also extended to a class of switching systems in continuous time. Based on these conditions, a systematic procedure for constructing discrete time switching sequences that, in a sense, globally asymptotically stabilise the system is then developed. The results are finally extended to systems with inputs in the sense of input-to-state stability and demonstrated on a small, benchmark example.

"Tradeoffs between quality-of-control and quality-of-service in large-scale

nonlinear networked control systems" by D.P. Borgers, R. Geiselhart and W.P.M.H. Heemels This paper investigates the properties of controlled dynamical systems where a communication network is included in the feedback loop and used to carry sensing and actuation information. The goal is to quantify the relation between the quality of service guarantees provided by the communication network and the performance of the closed loop system. In this paper quality of service is quantified through the maximum allowable transmission interval (a guarantee on the maximum time between two consecutive transmissions) and the maximum allowable delay (a similar guarantee on the maximum time between the transmission and reception of information). Control performance on the other hand is characterised through guarantees of input-to-state stability and the corresponding convergence rates; as usual a multiple Lyapunov function approach is adopted for this characterisation. The paper develops precise conditions on the effect of quality of service on stability, allowing one to determine bounds on the possible degradation of control performance if the quality of service deteriorates. The analysis is carried out in continuous time, starting with the very general framework of hybrid differential inclusions.

2.8. "Model measuring for discrete and hybrid systems" by T.A. Henzinger and J. Otop Within this article, the authors address the problem of measuring a model in the sense of determining the maximum permissible amount of perturbation applicable to constants occurring in the model that would guarantee preservation of satisfaction of its specification. This generalization of traditional Boolean verification is pursued for a broad set of models, ranging from weighted word and tree automata to parametric monotonic hybrid automata. As model measuring my be applied to systems not satisfying their specification, in which case model measuring would determine the minimum distance to a model variant satisfying the specification, the related question of model repair towards requirements satisfaction is also investigated.

2.9. "First steps toward formal controller synthesis for bipedal robots with experimental implementation" by A.D. Ames, P. Tabuada, A. Jones, W.-L. Ma, M. Rungger, B. Schrmann, S. Kolathaya and J.W. Grizzle In this paper, the authors apply formal methods to controller synthesis for bipedal robots. The control synthesis approach consists of two steps. Firstly, traditional nonlinear control techniques are used to stabilize the system on a low dimensional manifold; the evolution of the system on this manifold is described by the hybrid zero dynamics. Secondly, abstraction based control is applied to the hybrid zero dynamics to enforce safety specifications to be satisfied by the original system. The approach is validated using numerical simulations. Experimental results obtained by implementing the synthesized controllers on a bipedal robot are promising and show the relevance of the formal methods in control applications.

"A Type

  -based analysis of causality loops in hybrid systems modelers" by A.Benveniste, T. Bourke, B. Caillaud, B. Pagano and M. Pouzet This article addresses the static detection of causality loops in hybrid system models. It establishes a type system expressing dependencies between signals and enforcing absence of causality loops as a necessary condition for a system being well-typed. This ensures their causal consistency as necessary for avoiding deadlock situations during sequential simulation and for generating statically scheduled code.Martin Fränzle Department of Computing ScienceCarl von Ossietzky Universität Ammerlnder Heerstr. 114-118, 26129 Oldenburg, Germany Email adress: martin.fraenzle@informatik.uni-oldenburg.de Antoine Girard Laboratoire des Signaux et Systèmes CNRS, CentraleSupélec, Univ. Paris-Sud, Univ. Paris-Saclay 3 rue Joliot Curie, 91192 Gif-sur-Yvette, France Email adress: antoine.girard@l2s.centralesupelec.fr

	John Lygeros
	Automatic Control Laboratory
	ETH Zurich
	Physikstrasse 3
	CH-8092 Zurich, Switzerland
	Email adress: jlygeros@ethz.ch
	Sriram Sankaranarayanan
	Department of Computer Science
	University of Colorado, Boulder
	Boulder, CO, 80301, USA
	Email address: srirams@colorado.edu


Jonathan Siau
email: jonathan.siau@cedrat.com

Gérard Meunier
email: gerard.meunier@g2elab.grenoble-inp.fr

Olivier Chadebec
email: olivier.chadebec@g2elab.grenoble-inp.fr

Jean-Michel Guichon
email: jean-michel.guichon@g2elab.grenoble-inp.fr

Rémy Perrin-Bit
email: remy.perrinbit@cedrat.com

J.-M Guichon

Volume Integral Formulation Using Face Elements for Electromagnetic Problem Considering Conductors and Dielectrics

Keywords: Dielectric materials, magneto-harmonic, unstructured-PEEC, volume integral equation

A new integral formulation is presented, enabling the computation of resistive, inductive, and capacitive effects considering both conductors and dielectrics in the frequency domain. The considered application allows us to neglect any propagation effects and magnetic materials. In this paper, we will show how to improve the unstructured-partial element equivalent circuit approach to consider dielectric materials, keeping the same benefits. Results obtained with this formulation are compared to results from an industrial finite-element method software and measurements.

I. INTRODUCTION

F OR solving Maxwell's equations, many numerical methods and formulations can be found depending on the restraint assumptions. Finite-element methods (FEMs) are very popular for solving electromagnetic problems [START_REF] Albanese | Finite element methods for the solution of 3d eddy current problems[END_REF]- [START_REF] Dular | Dual finite element formulations for lumped reluctances coupling[END_REF], allowing us to assemble only sparse matrices, so any memory issues can be avoided. On the other hand, integral equation methods [START_REF] Albanese | Integral formulation for 3d eddy-current computation using edge elements[END_REF]- [START_REF] Rubinacci | Automatic treatment of multiply connected regions in integral formulations[END_REF] (IEMs) lead to assembling and storing full matrices. In the past, full matrices were really problematic due to the technological limitations. But the main advantage of IEM is that we need to mesh only the active regions, whereas the FEM requires meshing the surrounding air in addition. As a consequence, the IEM should be well suited to studying electromagnetic compatibility (EMC) problems. EMC applications mostly consist in studying distant components (such as wireless charging applications) or compact devices (such as printed circuit boards) and both have a large amount of air. So, using the FEM necessitates having a dense mesh in the air region around the devices, leading to practical issues (such as meshing difficulties or insufficient memory available). Also, the IEM enables us to compute the electromagnetic fields everywhere in space using only the active regions, so far interactions can be easily studied. Moreover, nowadays the IEM has regained interest due to new technologies and compression techniques. Many compression techniques exist, like clustering techniques such as the FMM [START_REF] Greengard | A fast algorithm for particle simulations[END_REF] or the H-matrices [START_REF] Hackbusch | A sparse matrix arithmetic based on H-matrices. I. Introduction to H-matrices[END_REF].

The applications of many electromagnetic formulations are often academic and their use for industrial purposes can be really a handful. A difficulty in modeling industrial applications in the electrical engineering domain is adding external electric components (like supply sources). A method has been proposed by Ruehli based on integral equations to solve efficiently electromagnetic problems coupled with an external circuit: the partial element equivalent circuit (PEEC) [START_REF] Ruehli | Equivalent circuit models for three-dimensional multiconductor systems[END_REF]. It consists in transforming a meshed electromagnetic device into its equivalent circuit made of lumped R-L-C elements and sources, so the coupling with the external circuit is natural. The classical PEEC was limited to structured meshes, but this limitation has recently been overcome by a Nonorthogonal PEEC formulation [START_REF] Ruehli | Nonorthogonal PEEC formulation for time-and frequency-domain EM and circuit modeling[END_REF] or by using face elements [START_REF] Whitney | Geometric Integration Theory[END_REF], [START_REF] Nguyen | An integral formulation for the computation of 3-d eddy current using facet elements[END_REF] for general meshes. This major improvement enables the treatment of more complex geometries and devices [START_REF] Bandinelli | A surface PEEC formulation for high-fidelity analysis of the current return networks in composite aircrafts[END_REF]. In practice, the PEEC method is mostly used to model conductors and resistive, inductive and capacitive effects [START_REF] Freschi | Unstructured PEEC formulation by dual discretization[END_REF]- [START_REF] Ferranti | Multipoint full-wave model order reduction for delayed PEEC models with large delays[END_REF]. A few previous publications have considered dielectric materials [START_REF] Heeb | Three-dimensional interconnect analysis using partial element equivalent circuits[END_REF], [START_REF] Ardon | MoM and PEEC method to reach a complete equivalent circuit of a static converter[END_REF], which have limitations. For instance, in [START_REF] Ardon | MoM and PEEC method to reach a complete equivalent circuit of a static converter[END_REF] the magneto-quasi-static and electrostatic problems are solved separately, so the phenomena are not strongly coupled. Moreover, the user has to define macroelements manually (which are a group of elements of the mesh) for the electrostatic solution, which may add some inaccuracy to the electrical phenomena.

The contribution of this paper is to propose a formulation which extends the PEEC method to take into account resistive, inductive, and capacitive effects in the presence of conductors and dielectric materials. Only one solution is needed here to solve the whole problem and the capacitive regions are supported by the mesh, so we have overcome the limitations of the previous literature. Moreover, the only assumptions here are that we neglect propagation effects (low frequencies) and that there is no magnetic material. Therefore, we will be using the freespace Green's function which does not take into account any delay. However, these assumptions are usual for low-frequency applications.

In the first section, we present the proposed formulation, starting from Maxwell's equations. Then, a practical formulation is introduced to overcome a limitation of the last one, in which a few assumptions will be stated. Last, the validation of the formulation is done with two examples. In the first one, we will compare the results of an industrial FEM software with the practical formulation, and in the second, the results are compared to real measurements. Note that the formulations will be expressed in the frequency domain, but they can be easily written in the time domain.

II. FORMULATION

Let us start by defining the following regions: Ω J and Ω D contain, respectively, the conductors and the dielectrics, and Ω = Ω J ∪ Ω D . Their borders are defined as Γ = ∂Ω, Γ J = ∂Ω J , and Γ D = ∂Ω D .

A. Maxwell's Equations

In a problem with conductors, dielectrics and without any magnetic material, we have the following constitutive laws:

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ J = σE in Ω J , D = E = 0 E + P in Ω D , B = μ 0 H (1)
where J is the current density, E the electric field, σ the conductivity, D the displacement current, P the polarization density, B the magnetic induction, the permittivity, and H the electromagnetic field. These quantities are expressed in SI units.

According to Maxwell's equations, we get

rotE = -jωB = -μ 0 jωH (2) rotH = J + jωD = (J + jωP) + 0 jωE (3)
with j the imaginary unit, ω = 2πf , and f the frequency. As in a vacuum, we can also write the potentials in the presence of dielectrics using the Lorentz gauge, as follows:

A = μ 0 4π Ω J + jωP r dΩ (4) jωV = 1 4π 0 Ω (J + jωP) ∇ 1 r dΩ (5)
where A and V are the magnetic vector potential and the scalar electric potential, respectively, and r is the distance between the testing and the integration points. We can note that only the jωP term has been added from the original expression.

A previous formulation was presented in [START_REF] Nguyen | An integral formulation for the computation of 3-d eddy current using facet elements[END_REF] to take into account resistive and inductive effects in the presence of conductors. The choice was made to use the current density J as the unknown. Moreover, as the current density's normal component is preserved between two elements, it seems natural to use the face element discretization. This also permits overcoming the structured mesh limitation of the classical PEEC. Now, we want to take into account capacitive effects and dielectric materials, so we have to use another unknown having the same conservative property in order to keep the same benefits. We choose the total current density J tot as the new unknown:

J tot = J + jωD = σ * E (6)
with σ * = σ + jω , which can be seen as a complex conductivity varying in the matter. J tot obviously keeps the same benefits in the conductors, and in the dielectric J tot = jωD verifies the same property. Thus, this integral equation formulation consists in matching the electrical constitutive law (6) everywhere in the materials (conductors and dielectrics) with E = -jωA -∇V [START_REF] Greengard | A fast algorithm for particle simulations[END_REF] resulting in

J tot σ * = -jωA -∇V. ( 8
)
To express the potential based on J tot , we use the second constitutive law and ( 6) to write

J tot = J + jωP + jω 0 E (9)
so we can deduce

J + jωP = σ * -jω 0 σ * J tot . ( 10
)
And by substituting ( 10) into (4), we can write (8) as

J tot σ * + jω μ 0 4π Ω σ * -jω 0 σ * J tot r dΩ = -∇V. ( 11
)

B. Discretization

To be able to treat unstructured meshes, we use the same test and trial functions as in [START_REF] Nguyen | An integral formulation for the computation of 3-d eddy current using facet elements[END_REF]. So, we discretize J tot using the face elements w j as follows:

J tot = k w k I k ( 12
)
with I j the current flowing through the face j.

Finally, applying a standard Galerkin projection to [START_REF] Whitney | Geometric Integration Theory[END_REF] with the same functions leads to the matrix system

([R t ] + jω[L t ]) {I} = {δV int } (13) [R t ] i,j = Ω w i w j σ * dΩ (14)
[L t ] i,j = jω Ω w i Ω σ * -jω 0 σ * w j r dΩdΩ (15)
(δV int ) i = - Ω w i ∇V dΩ (16)
with [R t ] a sparse matrix which represents the resistive effect in the conductor and capacitive effects in the dielectrics, and [L t ] is the inductive dense matrix. We can see an example of an equivalent circuit representation inFig. 2 of the mesh in Fig. 1. Note that Fig. 2 shows only the self-element contribution (not the mutuals), but the complete circuit can be easily constructed using [14, Fig. 1]. We are going to investigate the meaning of δV int , and we start by using the Divergence Theorem: Considering an inner face, we have

-(δV int ) i = Ω w i ∇V dΩ = Γ (w i • n)V dΓ - Ω (∇ • w i )V dΩ. (17)
(δV int ) i = Ω (∇ • w i )V dΩ (18)
and we know that the support of the face element w i is shared by two elements (e 1 and e 2 , see Fig. 1), so we have

(δV int ) i = e 1 V Vol(e 1 )
dΩ -

e 2 V Vol(e 2 ) dΩ (19)
knowing that ∇ • w i = Vol(e i ) -1 , with Vol(e i ) giving the volume of the element e i . We can see that [START_REF] Ardon | MoM and PEEC method to reach a complete equivalent circuit of a static converter[END_REF] represents the difference of averaged potential between two inner elements. So, from a circuit approach, we can assimilate this as the applied voltages. Let's consider a face j on the border. Its support is only on one element (e 1 for example), so [START_REF] Antonini | An accurate interpolation strategy for fast frequency sweep of partial element equivalent circuit models[END_REF] gives

(δV int ) j = Γ (w j • n)V dΓ - e 1 V Vol(e 1 ) dΩ (20)
and knowing that (w j • n) = S -1 j , with S j the surface of the face j, we have

(δV int ) j = Γ j V S j dΓ - e 1 V Vol(e 1 ) dΩ (21)
also representing the potential difference between a border face and its inner support.

So, the right-hand side vector is the potential difference between the dual mesh branches (see Fig. 1), where the dual mesh represents the equivalent electrical circuit.

Finally, (13) corresponds to the incomplete circuit equations (U = ZI): the circuit is not closed (see Fig. 1). So, we will have to add new branches on the border faces, and to express their contributions. For example, we have to add four capacitive branches in Fig. 2 (linking the nodes N 2 , N 3 , N 5 , and N 6 to the common node ∞) in order to close the circuit. To complete the equations, we have to express the difference of potential between the border faces, so we introduce the common node ∞ and

(δV ext ) j = V Γ j -V ∞ ( 22
)
with V Γ j and V ∞ the potential of the border face j and of the common node ∞, respectively. Actually, we have already expressed V Γ j in the previous case:

V Γ j = Γ j V S j dΓ j . ( 23
)
We suppose that the node ∞ has null potential (V ∞ = 0), which leads to

(δV ext ) j = V Γ j . (24)
Last, we can use [START_REF] Le-Van | A volume integral formulation based on facet elements for nonlinear magnetostatic problems[END_REF] to get

(δV ext ) j = 1 4π 0 Γ 1 S j Ω σ * -jω 0 jωσ * J tot ∇ 1 r dΩdΓ (25)
and using the same discretization procedure, we can write

{δV ext } = 1 jω [P t ]{I} (26)
with

[P t ] i,j = 1 4π 0 Γ 1 S i Ω σ * -jω 0 σ * w j ∇ 1 r dΩdΓ. (27)
Here, we have managed to express the difference of potential on the border faces, which can be assimilated as capacitive effects. Last, we have to solve the equation:

[R t ] + jω[L t ] + 1 jω [P t ] {I} = {δV } ( 28
)
with δV = δV int + δV ext the difference of potential of all the branches. We saw that (28) actually corresponds to a circuit equation, so we can solve it using a circuit solver. The main advantage is that this permits ensuring the conservation of current (Kirchhoff's first law) in the problem. The solver circuit we chose uses the independent loop search approach, instead of the modified nodal approach (usually picked in the literature [START_REF] Freschi | Fast block-solution of PEEC equations[END_REF]). This choice was made to minimize of the number of unknowns. However, the final equation to solve is

[M] [R t ] + jω[L t ] + 1 jω [P t ] [M] T I m = 0 (29)
with [M] the change of basis matrix from the branches to the independent-loop basis and I m the loop current vector.

At this point, we have presented a general approach to include dielectric materials and the capacitive effects. This formulation needs matrices with a dependency on the frequency (σ * and σ *j ω 0 σ * for instance). This necessitates assembling each matrix once for each frequency we want to study. So, for a frequency study, this is a critical disadvantage, and necessitates using a more practical formulation.

III. PRACTICAL FORMULATION

In this section, we will present a more practical formulation using the same approach, but making a few more assumptions to overcome the last critical point.

First, we assume that σ ω 0 (30) so the previous unknown J tot becomes

J tot = J in Ω J jωD in Ω D .
Following the same protocol as in Section II, we write the Ampère-Maxwell equation ( 7) as

E + jω μ 0 4π Ω J + jωP r dΩ = -∇V (31)
and we can also write

J + jωP = ⎧ ⎪ ⎨ ⎪ ⎩ J = J tot in Ω J jωP = -0 J tot in Ω D (32)
using [START_REF] Albanese | Finite element methods for the solution of 3d eddy current problems[END_REF]. In consequence, (31) can be rewritten as

E + jω μ 0 4π Ω J J tot r dΩ + Ω D -0 J tot r dΩ = -∇V.
(33) Now, we can apply the Galerkin projection with the interpolation functions w i :

Ω w i EdΩ + jω μ 0 4π Ω w i Ω J J tot r dΩ + Ω D -0 J tot r dΩ dΩ = - Ω w i • ∇V dΩ. ( 34
)
And we can easily see from the first two constitutive laws that

E = ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ J tot σ in Ω J J tot jω in Ω D
and, by using the same discretization as [START_REF] Nguyen | An integral formulation for the computation of 3-d eddy current using facet elements[END_REF], we finally have

[ R] + jω[ L] + 1 jω [ C] {I} = {δV int } (35) with [ R] i,j = Ω J w i w j σ dΩ (36) [ C] i,j = Ω D w i w j dΩ (37) [ L] i,j = μ 0 4π Ω w i Ω J w j r dΩ+ Ω D -0 w j r dΩ dΩ. ( 38
)
We can see that we have one more matrix than for (13) ([ C]). In fact, the previous matrix [R t ] has been simplified to

[ R] + 1 j ω [ C].
So, even if we have added a matrix, the actual amount of memory needed has been divided by 2, because [R t ] is complex while [ R] and [ C] are real matrices. We can make the same remark about [ L] and [L t ], so the whole storage has been globally divided by 2.

At this point, we still have to complete the circuit equation to close the equivalent circuit, and to take into account the capacitive effects. The protocol is exactly the same, so let's begin with [START_REF] Duffy | Feature selective validation (FSV) for validation of computational electromagnetics (CEM). Part I-The FSV method[END_REF]:

(δV ext ) j = V Γ j = Γ j V S j dΓ j . ( 39
)
Using (32) in ( 5) leads to

(δV ext ) j = 1 4π 0 1 jω Γ 1 S j Ω J J tot ∇ 1 r dΩ + Ω D -0 J tot ∇ 1 r dΩ dΓ. (40)
In this case, we can simplify the first twofold integral:

Γ 1 S i Ω J J tot ∇ 1 r dΩdΓ = Γ 1 S i Γ J ( J tot • n) 1 r dΓdΓ (41)
using the Divergence Theorem and knowing that div J tot = 0. In addition, assuming that the dielectrics are linear (which implies divP = 0 in Ω D ), we can reapply the Divergence Theorem to simplify the second twofold integral

Γ 1 S j Ω D -0 J tot ∇ 1 r dΩdΓ = Γ 1 S j Γ D -0 ( J tot • n) 1 r dΓdΓ. ( 42
)
Now, we use the same discretization on J tot and we can write

{δV ext } = 1 jω [ P]{I} (43)
with

[ P] i,j = 1 4π 0 Γ 1 S i Γ J 1 S j 1 r dΓ + Γ D -0 1 S j 1 r dΓ dΓ. ( 44
)
Here, we have managed to complete the circuit equations with a matrix considering only boundaries. This matrix is much smaller than [P t ], so it requires less memory and can be assembled faster. At the end, we have the next equation

[ R] + jω[ L] + 1 jω [ C] + [ P] {I} = {δV } (45)
to solve. We can note that the matrices do not have any dependency on the frequency, so we have overcome the last critical point. Moreover, we managed to save more than the half of the memory required before. And also, we note that the matrices are purely real, so, using real-assembly codes permits lowering the assembly time.

We also note that the full matrices ([ L] and [ P]) are classical in the theory of integral equations and can be compressed. So, this could allow us to solve bigger problems.

IV. VALIDATION

The validation of the practical formulation is done comparing the impedances obtained for two examples: a microcoil and a real device. The first example will compare our results to those from an industrial software, and the second will compare our results against measurements.

A. Case 1: Microcoil

The considered microcoil (see Fig. 3) is made of three conductive spires with a dielectric between the I/O. This validation is done comparing the impedances from the formulation and an industrial finite element software: Flux 3D [START_REF]Flux R 3D Software[END_REF]. Flux3D takes advantage of the symmetry plane on XY , while our integral code cannot. So, the results will show only a half of the device on Figs. 4 and5, even though the practical formulation used the entire mesh. For a matter of ease, the Flux3D uses a structured mesh for the active regions (conductor and dielectric) composed of 9579 hexahedra and an usual mesh is used in the air with 92 049 elements. For the proposed formulation, we used a mesh of 4822 tetrahedra.

The solution of this formulation was done using an FGMRes solver and a Block lower-upper (LU) preconditioner. The FEM solutions used an ICCG solver, so the computational time cannot be fairly compared. The computations have been made using a Dell Precision M4800 with an Intel(R) Core(TM) i7-4800MQ CPU @2.70 GHz and 32Go RAM.

The numerical results are presented in Fig. 6 for the impedance modulus and in Fig. 7 for the impedance phase. In the last two figures, "FEM" stands for the results provided by Flux 3D, which are published in [START_REF] Pham Quang | Validation of a 3d decoupled magnetodynamic-electric model by determining the impedance of microcoils[END_REF].

According to the IEEE standard [START_REF]Standard for Validation of Computational Electromagnetics Computer Modeling and Simulation-Part 1[END_REF], the comparison and the qualification of the last curves is done using the FSV techniques [START_REF] Duffy | Feature selective validation (FSV) for validation of computational electromagnetics (CEM). Part I-The FSV method[END_REF], [START_REF] Orlandi | Feature selective validation (FSV) for validation of computational electromagnetics (CEM). Part II-Assessment of FSV performance[END_REF] via the stand-alone software FSV tool [START_REF] Orlandi | Free stand-alone FSV application[END_REF]. We applied the amplitude difference measure (ADM) and the feature difference measure (FDM) techniques to the impedance modulus curves [see Fig. 8 We can see in Figs. 6 and 7 that the results of both methods are very similar. This is also supported by the FSV technique

B. Case 2: Real Device

Let's now consider a real device made of two layers of per 35 μm thick separated by a layer of dielectric FR4 with a thickness of 1.47 mm (see Fig. 10). Two cases have been considered to do the measurements: The short-circuit case consisting of connecting each layer by a corner (see Fig. 11) and the open-circuit case which consists in leaving the layers without a connexion between them (see Fig. 12). Both cases are considered from 1 to 100 MHz.

The same mesh has been used for all the studies and is made of 3210 hexahedra which leads to 16 951 faces/degrees of freedom. We are using a structured mesh because the geometry is wellsuited and the conductors are very thick, but an unstructured mesh could also have been used. The thickness of each layer is discretized with only one element, because the skin depth is greater than the thickness for the studied frequencies.

The solutions were done using a classical LU factorization, due to the ill-conditioning of the problem. The computations were done on the same laptop as in the previous section. The numerical results are presented in Figs. 13 and 14 for the shortcircuit and the open-circuit case, respectively.

On one hand, both cases show a good accuracy regarding the first peak frequency. On the other hand, the amplitude of the peaks are underestimated. But overall, we see a good match between the curves. To quantify this, we have used the same FSV technique as before (using the FSV Tool). The results are shown in Fig. 15

V. CONCLUSION

We have presented a new formulation to consider resistive, inductive, and capacitive effects in the presence of conductors and dielectrics, with very few assumptions. We pointed out that its generality also spoils its efficiency, due to the frequency dependency of its matrices. So, to solve a problem for multiple frequency one should use the second formulation.

To overcome this dependency, we have presented the practical formulation with a few more assumptions. We have pointed out that these additional restrictions allow us to assemble simpler matrices and to greatly reduce the required storage compared to the general formulation.

The practical formulation has been validated considering two cases. The first consisted in comparing the impedance of a virtual device against an industrial FEM software. The results of both methods show a great match and have been qualified as excellent by the software FSV Tool. The last case consisted in comparing the computed impedances against measurements on a real devices, and the curves show a great match, rated as good to excellent. Further research will compare the proposed formulation to others regarding the computational cost and time.

Fig. 1 .

 1 Fig.1. Example of a topology at the border between a conductor and a dielectric. Each face of the primal mesh corresponds to a dual branch and the black (circuit) nodes correspond to volume elements and border faces.

Fig. 2 .

 2 Fig. 2. Simplified equivalent circuit (the mutuals are hidden) of the geometry shown in Fig. 1.

Fig. 3 .Fig. 4 .

 34 Fig. 3. Studied device: Microcoil with three spires.

Fig. 5 .Fig. 6 .

 56 Fig. 5. | J tot | isovalues at 100 MHz using the "practical formulation."

  (a) and (b)] and to the impedance phase curves [see Fig. 9(a) and (b)].

Fig. 7 .

 7 Fig. 7. Impedance phase versus frequency. In blue (dashed): The reference (FEM), in red: the simulation results.

Fig. 8 .

 8 Fig. 8. (a) Results of the ADM applied to the impedance modulus curves. (b) Results of the FDM applied to the impedance modulus curves.

Fig. 9 .

 9 Fig. 9. (a) Results of the ADM applied to the impedance phase curves. (b) Results of the FDM applied to the impedance phase curves.

Fig. 10 .

 10 Fig. 10. Picture and dimensions of the studied device.

Fig. 11 .

 11 Fig. 11. Current density isovalues of the short-circuit case at 1 MHz, with the circuit connexions drawn. The modulus | J tot | and the real part ( J tot ) of the current density are represented by the colored elements and vectors, respectively.

Fig. 12 .

 12 Fig. 12. Current density isovalues of the open-circuit case at 1 MHz, with the circuit connexions drawn.

  (a) and (b) for the open-circuit case and in Fig.16(a) and (b) for the short-circuit case. The ADM technique rates both case mostly from good to excellent, whereas the FDM rates them mostly as excellent.

Fig. 13 .

 13 Fig.[START_REF] Bandinelli | A surface PEEC formulation for high-fidelity analysis of the current return networks in composite aircrafts[END_REF]. Impedance modulus versus frequency for the short-circuit case. In blue (dashed): The reference (measurements), in red: the simulation results.

Fig. 14 .

 14 Fig. 14. Impedance modulus versus frequency for the open-circuit case. In blue (dashed): The reference (measurements), in red: the simulation results.

Fig. 15 .

 15 Fig. 15. (a) Results of the ADM applied to the impedance modulus curves for the open-circuit case. (b) Results of the FDM applied to the impedance modulus curves for the open-circuit case.

Fig. 16 .

 16 Fig. 16. (a) Results of the ADM applied to the impedance modulus curves for the short-circuit case. (b) Results of the FDM applied to the impedance modulus curves for the short-circuit case.

Jonathan SIAU was born in 1989. He received the Master's degree in computational science and applications from the University of Rennes 1, Rennes, France, in 2013. Since 2014, he has been working toward the Ph.D. degree entitled New solvers for 3-D power electronics devices simulation, thanks to a collaboration of CEDRAT S.A and the Grenoble Electrical Engineering Laboratory.

His research interests include low-frequency formulations development for the power electronics, their numerical properties, and numerical methods. He is currently the CNRS Senior Researcher (directeur de recherche CNRS) and leads the Models, Methods, and Methodologies Applied to Electrical Engineering Research Group, G2ELab, University of Grenoble, Grenoble, France. He has been a Visiting Scientist at the Technology Centre of the Federal University of Santa Catarina, Brasil, during 2012-2013. His research interests include computational electromagnetism (finite element and integral methods), equivalent magnetic sources identification by solving inverse problems, and low magnetic fields measurements. He has coauthored more than 140 papers published in international journals and conference proceedings. He is currently System, Physics, and Application Manager of Cedrat, Meylan, France. His activities focus on computational electromagnetism (finite element) in Flux Cedrat software.

Gérard

Jean-Michel


Arnaud Bourdin

Anne Sophie Gamez

Isabelle Vachier

Bruno Crestani

LAM is another small airway disease: lessons from microCT