Ezziane Mohammed

Kadri Tahar

Molez Laurent
email: laurent.molez@insa-rennes.fr

Jauberthie Raoul

Belhacen Ali

Laurent Molez

High temperature behaviour of polypropylene fibres reinforced mortars

Keywords: Polypropylene fibre mortar, Mechanical behaviour, Fracture energy, Fire test, Spalling, Punching shear test

The aim of this paper is primarily experimental and is intended to analyse the behaviour of two cementitious materials, before and after heat treatment: one unreinforced (i.e. without fibres) and the other reinforced (with polypropylene fibres).

At room temperature and after heating up to 500°C, the bending strength is improved by the presences of fibres. The residual young modulus is slightly higher for the fibres reinforced samples.

As the temperature increases, the strength gain due to fibres inclusion is reduced. Beyond 500°C, the bending strength is lower for the fibre reinforced cementitious material compared to those without fibres.

Fracture energy is also improved for the fibre mortars at room temperature. At 400°C this improvement decreases gradually with the introduction of polypropylene fibres. Beyond this temperature and due to the introduction of polypropylene fibres, the fracture energy is reduced.

Another test is developed: rapid heating due to exposure to a flame. The temperature in the front side reaches in few seconds 1000°C. At this temperature and after one hour of exposure, the opposite side reached 140°C. After cooling, the punching shear strength of the fibre mortar is definitely weaker than of the mortar without fibre.

INTRODUCTION

The polypropylene fibres are used in the mortars and concrete to decrease the plastic shrinkage, cracking and micro-cracking of surface [START_REF] Grzybowski | Shrinkage cracking of fiber reinforced concrete[END_REF]. They increase cohesion and reduce slump [START_REF] Johnston | Proportioning, mixing and placement of fibre-reinforced cements and concretes, Production Methods and Workability of Concrete[END_REF]. While inclusion of metal fibres improves the fire resistance [START_REF] Lie | Thermal and mechanical properties of previous steel fibre-reinforced concrete at elevated temperature[END_REF], the presence of polypropylene fibres in the material does not [START_REF] Serdar | High temperature resistance of normal strength and autoclaved high strength mortars incorporated polypropylene and steel fibres[END_REF].

Only the effect limiting spalling at high temperatures is about to be unanimously recognized. This action is primarily studied in high performance concrete i.e. with low water/cement ratio [START_REF] Serdar | High temperature resistance of normal strength and autoclaved high strength mortars incorporated polypropylene and steel fibres[END_REF][START_REF] Bilodeau | Optimization of the type and amount of polypropylene fibres for preventing the spalling of lightweight concrete subjected to hydrocarbon fire[END_REF][START_REF] Nishida | Study on the properties of high strength concrete with short polypropylene fibre for spalling resistance, Concrete under severe conditions: environment and loading[END_REF].

However, during a fire, very strong thermal gradients take place in first centimetres of concrete and the thermal damage rapidly decreases from a maximum to nil [START_REF] Colombo | New NDT for techniques for the assessment of fire-damaged concrete structures[END_REF][START_REF] Mindeguia | Temperature, pore pressure and mass variation of concrete subjected to high temperature -Experimental and numerical discussion on spalling risk[END_REF]. So, in these first centimetres, where temperatures are less than a critical temperature, polypropylene fibres should still have an effect on mechanical behaviour. The objective of this work is then to examine the effect of the temperature on the residual behaviour of the polypropylene fibre mortars. In addition to the classical flexural strength, we focus to the cracking behaviour and more particularly to the fracture energy and to the stress intensity factor. These characteristics are important in connection with the crack resistance and in particular with the resistance to spalling.

In this paper two type of heating were studied: low rate heating (2°C/min) and high rate heating (rate similar to ISO 834 standard fire). These two studies allow to separate physicochemical effects and temperature gradient effect.

MATERIALS

To ensure an effective bridging of the cracks (to limit the crack opening), the dimension of the selected fibres must be greater than those of the aggregates [START_REF] Serdar | High temperature resistance of normal strength and autoclaved high strength mortars incorporated polypropylene and steel fibres[END_REF]. At room temperature, the aim of the polypropylene fibres in the composite matrix is to ensure a macro cracking bridging and to maintain high post peak strength at a very large crack [START_REF] Hannant | Fibre cement and fibre concrete[END_REF]. In this study cementitious mortars were studied, and then, 12 mm monofilament polypropylene fibres were chosen. Characteristics of the fibres (given supplier: SIKA) are detailed in table 1. The mortars are made with standardized sand according to CEN 196-1 standard with the ISO 679. The mass proportions of cement, sand and water are 1:3:0.5. The volume proportion of added fibres in the fibre mortars is 0.58% (i.e. 5.2 kg/m 3 ). Water was added to cement and mixed to obtain a homogeneous paste; sand then gradually added to the paste and mixed until homogeneous. Fibres were added at the final stage and dispersed manually. The constituents were mixed for two minutes after the introduction of the fibres.

The specimens were then stored in a wet room (20°C, 95% RH) for 7 days and then stored in dry room (20°C, 50% RH) up to an age of 28 days. Under these conditions, a significant portion of free water in the cement matrix had evaporated [START_REF] Andrade | Relative humidity in the interior of concrete exposed to natural and artificial weathering[END_REF].

SLOW HEATING

Heat exposure in an oven

Samples, 4x4x16 cm 3 , were heated in an electric furnace to the desired temperature at a rate of 2°C/min.

The exposure temperature was maintained during one hour (1h) and cooling to room temperature was carried out in the closed and disconnected furnace (approximately -0.3°C/min). The controlled temperature is measured in oven chamber (i.e. not in samples). This process conducts to low thermal gradients. The damages in mortar are mainly of physicochemical origin [START_REF] Georgali | Microstructure of fire-damaged concrete, A case study[END_REF].

Flexural strength

After cooling, samples were then tested in a four-point bending configuration. Six specimens were tested for each condition (mortar type and heating temperature). The load is measure with a 50 kN load cell. A template is attached to the specimen to measure deflection using two LVDTs (+/-1mm) (figure 1). The bending test results carried out on unreinforced mortar (MN) and fibre reinforced mortar (MNP) are presented in figure 2. The post peak residual load clearly highlights the role of fibres when samples are not heated. The stress transfer between the faces of the cracks is significant and result is an increased ductility. As already observed by other authors [START_REF] Pliya | Influence des fibres de polypropylène sur le comportement de bétons soumis à une température élevée. Colloque International sur la caractérisation et modélisation des matériaux et structures[END_REF][START_REF] Xiao | On residual strength of high-performance concrete with and without polypropylene fibres at elevated temperatures[END_REF][START_REF] Behnood | Comparison of compressive and splitting tensile strength of highstrength concrete with and without polypropylene fibres heated to high temperatures[END_REF], this phenomenon attenuates very quickly between 400 and 500°C and disappear beyond 500°C. Figure 3 represents the evolution of the bending strength with respect to the exposure temperature. The polypropylene fibres have a positive role until 400°C as it was stated previously. The experimental values are reported in table 3.

Young modulus

Young modulus is defined as the elastic (and linear) stage of the load-deflection curve. The evolution of the residual Young modulus with respect to the exposure temperature is shown in figure 4 and table 4.

According to the test results, we notice that the modulus of the sample with fibre is greater than that of the non-fibre samples. The variation is accentuated until 400°C. At 500°C; they are identical and beyond this temperature; the young modulus of the fibred samples becomes almost zero.

Fracture energy

The fracture energy is defined as the area under the stress-strain curve of the four-point bending test [START_REF] Guo | Determination of quasibrittle fracture law for cohesive crack models[END_REF].

Figure 5 and table 5 show the test results of the fracture energy with respect to the exposure temperature.

When the samples are not heated, the fracture energy of the fibre specimen is about three times the energy of the non-fibre sample. This difference in fracture energy tends to converge and nullify at about 500°C.

Above 500°C, gain increase in fracture energy decreases for samples with fibres.

Stress intensity factor

The stress intensity factor K I characterizes the resistance of material to the propagation of the crack and to the damage [START_REF] Charmet | Mécanique du solide et des matériaux: Elasticité-Plasticité-Rupture[END_REF][START_REF] Menou | Residual fracture energy of cement paste, mortar and concrete subject to high temperature[END_REF][START_REF] Pierre | A finite element study of the influence of micro-fibres on the tensile strength of cement paste[END_REF]. This parameter can be deducted from the calculation of the fracture energy G f and of the Young modulus E. It is defined by the following relation:

  0.5 If K G .E  (eq. 1)
This parameter takes into account the degradation of the cement matrix; rather quantifiable by the young modulus and the fibre degradation (modification of the characteristics of fibres, loss of cohesion), rather quantifiable by the fracture energy. This is thus a representative parameter of the overall damage of the fibrous material.

The evolutions of the stress intensity factor of the mortars with respect to the exposure temperature are illustrated in figure 6 and table 6. According to the test results (figure 5 and figure 6), we notice that the fracture energy and the stress intensity factor have similar behaviour when exposed to equivalent temperatures.

The stress intensity factor for the fibre mortar K I has a great value at 20°C and 400°C. Beyond this limit, we observe a very pronounced decreasing of the stress intensity factor. This decrease is mainly due to the fact that the polypropylene fibres have certainly totally melted and calcined. The fibres that stop crack propagation are no longer present.

At the high temperatures, two types of damage appear: (1) the physicochemical modifications and (2) cracking of the cementing matrix [START_REF] Georgali | Microstructure of fire-damaged concrete, A case study[END_REF][START_REF] Fu | Experimental study of micro/macro crack development and stress-strain relations of cement-based composite materials at elevated temperatures[END_REF][START_REF] Fares | Self-consolidating concrete subjected to high temperature[END_REF]. Results from a previous study showed that the polypropylene fibres had more effect on the reduction of cracking than on the reduction of the physicochemical modifications [START_REF] Ezziane | Heat exposure tests on various types of fibre mortar[END_REF]. From 170°C, the polypropylene fibres melt and create a connected porosity [START_REF] Komonen | Effect of high temperature on the pore structure and strength of plain and polypropylene fibre reinforced cement pastes[END_REF]. This porosity allows evacuating vapour over-pressures and thus reducing cracking [START_REF] Kalifa | High-temperature behaviour of HPC with polypropylene fibres from spalling[END_REF].

This weaker cracking explains the best results in term of resistance, ductility and fracture energy or stress intensity factor.

Moreover, when the temperature of the furnace is 400°C, the temperature inside the specimen is lower and can be not sufficient to destroy all fibres. The propagation of heat is slowed down in the dehydrated cementing matrix [START_REF] Georgali | Microstructure of fire-damaged concrete, A case study[END_REF][START_REF] Pliya | Influence des fibres de polypropylène sur le comportement de bétons soumis à une température élevée. Colloque International sur la caractérisation et modélisation des matériaux et structures[END_REF][START_REF] Xiao | On residual strength of high-performance concrete with and without polypropylene fibres at elevated temperatures[END_REF][START_REF] Behnood | Comparison of compressive and splitting tensile strength of highstrength concrete with and without polypropylene fibres heated to high temperatures[END_REF][START_REF] Guo | Determination of quasibrittle fracture law for cohesive crack models[END_REF][START_REF] Charmet | Mécanique du solide et des matériaux: Elasticité-Plasticité-Rupture[END_REF][START_REF] Menou | Residual fracture energy of cement paste, mortar and concrete subject to high temperature[END_REF]. The fibres disappear gradually from the periphery of the sample towards the centre: the kinetics of propagation of heat is slowed down by the more or less insulating dehydrated cementing matrix [START_REF] Mindeguia | Temperature, pore pressure and mass variation of concrete subjected to high temperature -Experimental and numerical discussion on spalling risk[END_REF]. Fibres remain in the heart of the sample and contribute to the ductility of material.

Beyond 500°C, the fibres completely disappeared to create an additional porosity [START_REF] Komonen | Effect of high temperature on the pore structure and strength of plain and polypropylene fibre reinforced cement pastes[END_REF][START_REF] Kalifa | High-temperature behaviour of HPC with polypropylene fibres from spalling[END_REF][START_REF] Noumowé | Transit heating effect on high strength concrete[END_REF][START_REF] Hager | Residual transport properties of heated HPC modified with polypropylene fibres, 1st International Workshop on Concrete Spalling due to Fire Exposure[END_REF][START_REF] Min | Mechanical properties of high-strength concrete after fire[END_REF]. This porosity leads thus to a lower energy than that of the not fibred samples. This process has been observed by optical microscopy (figure 7). Small slice extracted from different depth of a sample exhibits completely calcined polypropylene fibres in the first millimetres. Deeper slice contains yet partially melted fibre.

FLAME TEST

Exposure to fire

The second heat treatment consists of placing plates of dimensions 16x16x4 cm 3 made of unreinforced mortar and polypropylene fibre reinforced mortar in a flame test apparatus (figure 8a). Five thermocouples were inserted in the plate (figure 8b) at different depth and two other thermocouples were positioned in contact with the internal face (hot face) and external face (cold face) (figure 8c). The temperature of the exposed face (dimensions 16x16 cm 2 ) is forced to 1000°C and was maintained for one hour. Figure 9 compares the rate of temperature between the experimental curve and the ISO 834 standard fire. The experimental rate of temperature is similar to that of ISO 834 up to 1000°C. This rate leads to a strong temperature gradient in the sample. The sample dimensions are assumed to be enough to observe spalling.

The heat flux is assumed to be unidirectional. The flame position allows setting in a precise manner the temperature of the surface.

(a) (b) (c) Two series of test were performed:

-First series (equipped with 7 thermocouples (type K)): the evolution of the temperature of this series is obtained by placing thermocouples (drilled after casting) at different depths of the sample (external face -1 cm -1.5 cm -2 cm -2.5 cm -3 cm, and the internal side) as illustrated in figure 8c.

-Second series (equipped with 2 thermocouples): thermocouples were positioned in contact with the external (cold face) side and at the internal (hot face) side. It has to be noted that this thermocouple measures radiation and convection from the flame, rather than surface temperature at the exposed face.

However it was very difficult to achieve the measuring of the surface real temperature.

At the end of the test, cooling was performed naturally at room temperature. Punching shear tests were, then carried out for a comparison.

Heat transfer

During these tests, and after 5 minutes of treatment, we observed, drops of water and water vapour escaping from the drilling of thermocouples of the unexposed surface of the test sample for the 2 types of mortars. Nine minutes after treatment, we observed the apparition of a fist crack in the standard mortar.

This crack has developed from lateral side to the centre. However, this phenomenon did not appear for fibre reinforced mortar.

Figure 10 shows the evolution of the temperature of both mortars with respect of time and depth. The temperature evolution consists of three stages: a rapid increase of the temperature until 100°C, an isothermal stage around 100°C and an increase of the temperature. The isothermal stage around 100°C corresponds to the evaporation of free water at different depths. Mindeguia et al., [START_REF] Mindeguia | Temperature, pore pressure and mass variation of concrete subjected to high temperature -Experimental and numerical discussion on spalling risk[END_REF] observed this plateau at 100°C for low compactness concrete.

Mass loss after tests

The mass loss at the end of the test is almost identical for the two mortars. The mass loss of the standard mortar (MN) and fibre mortar (MNP) is 6.4% and 6.5% respectively. We can say that the fibres embedded in the mortar have virtually no influence on the evaporation of water from mortar and hydroxides OH after this heating. Similar observations were reported by [START_REF] Pliya | A way to improve the behaviour of concrete at high temperature: addition of a cocktail of polypropylene and steel fibres[END_REF].

Residual strength

Subsequently, to enable a comparison between the residual mechanical strength of samples a punching shear strength [START_REF] Menétrey | Synthesis of punching failure in reinforced concrete[END_REF] test was performed on non-drilled samples as presented in figure 12. Results are reported on table 7. Each value is the mean between ten samples; the difference between these ten values is maximum 3%. Before flame exposure, punching shear strengths are quite similar for the two mortars.

After flame exposure, all mortars show a loss of punching shear strength. Fibre mortars have a lower residual strength compared to those without fibres: strength reduction is 64.5% for mortar without fibre and 73.7% for mortar with polypropylene fibres. A temperature of 1000°C at the hot face for 1 hour conducts to temperature greater than 170°C (polypropylene melting) up to about 1 cm from the cold face:

it can be considered that all fibres were completely melted up to about 1 cm from the cold face. The induced porosity due to the degradation of the fibres decreases the strength. However, after heating, a part of polypropylene fibres is still present in the outer face side. The cold face temperature is only 140°C

after one hour, while the melting point is given at 170°C. Thus, the strength loss is limited. These results are consistent with the results obtained in a preliminary study [START_REF] Ezziane | Evolution of the punching strength of mortar subjected to fire[END_REF].

Spalling during heating

Induced porosity due to melting of polypropylene fibres is considered to facilitate the transfer of steam and avoid spalling [START_REF] Kalifa | High-temperature behaviour of HPC with polypropylene fibres from spalling[END_REF][START_REF]Eurocode, Design of concrete structures[END_REF][START_REF] Khoury | Polypropylene fibres in heated concrete. Part 1: molecular structure and materials behaviour[END_REF]. It is reasonable to consider that concrete incorporating polypropylene fibre can provide a benefit to concrete so as to prevent it from explosive spalling, due to the fact that it is melted under temperature around 170°C and hence accumulated moisture pressure in concrete can escape through inter-connected pores to outside of concrete. Samples from mortars with and without fibre at a water to cement ratio equal to 0.5 did not spall. However, other tests with polypropylene fibre reinforced mortars with the same fibre volume ratio but a water/cement ratio of 0.4 spalled during testing (figure 13).

Similar observations were conducted in studies on concrete with low water/cement ratio (concrete with 0.26 water/cement ratio and 1 kg/m 3 of polypropylene fibres spalled, but concrete with 0.26 water/cement ratio and 0.6 kg/m 3 of polypropylene fibres and 40 kg/m 3 of steel fibres did not spall [START_REF] Peng | Explosive spalling and residual mechanical properties of fibre-toughened high-performance concrete subjected to high temperatures[END_REF]) or on lightweight concrete (concrete with 0.33 ratio spalled, and concrete with 0.42 ratio did not spall [START_REF] Bilodeau | Optimization of the type and amount of polypropylene fibres for preventing the spalling of lightweight concrete subjected to hydrocarbon fire[END_REF]).

CONCLUSION

This study highlights the contribution of polypropylene fibres in cement based building materials exposed to an increasing temperature or exposed to flame. We have studied the mechanical behaviour of standard mortars and fibre reinforced mortars heat-treated at 400°C, 500°C, 600°C, 700°C, 800°C and 1000°C.

Two types of treatment have been made: a heating furnace (low rate) and an extreme heated test with flame (high rate).

MNP

At room temperature (before heating), the post peak stress clearly highlights the role of the fibres: the stress transfer between the faces of the cracks by polypropylene fibres is important and gives samples ductility.

After heating, below 400°C, the fibres have an uncracking effect by allowing the dissipation of fluid over pressure in the matrix. This phenomenon disappears at a temperature over 500°C. The fibre samples have greater young modulus compared to that of the non-fibred ones with an increasing gap until a temperature of 400°C and decreases beyond this value. At 500°C, they are identical. After 500°C, the young modulus becomes almost zero. The fracture energy at 20°C is 3 times greater for the fibred samples compared to the non-fibred ones. This ratio drops to 1.5 at 400°C. Fibres effects are negative beyond the temperature of 500°C.

It appears, from the flame test that the fibre mortar and non-fibre mortar have the same thermal behaviour. After one hour, the unexposed side is about 140°C.

For the residual punching strength; the results show that polypropylene fibres are not so efficient: the polypropylene fibres limit the cracking during heating but do not increase the punching strength.

By melting at 200°C the polypropylene fibres are considered to create a porosity which allows a limitation of the pressure due to the evaporation of water treatment and thus a limitation of cracking. The spalling phenomenon was not observed for all mortars having a water/cement ratio of 0.5. However, spalling was observed for mortars with water/cement ratio of 0.4.

For design of concrete structures, it has to be noted that the temperature remains below 500°C after two centimetres inside the concrete (i.e. after the cover over reinforcing steel). Up to 500°C, polypropylene reinforced mortar gives a pseudo-ductile behaviour with a better strength, a better Young modulus, a better fracture energy than mortar without fibre. Thus, the safety of concrete structures submitted to fire is increased. After fire, the rehabilitation of concrete cover will be limited and associated cost decreased.

Figure 1 :

 1 Figure 1: Four-point sample with measuring t.

Figure 2 :

 2 Figure 2: Representative curves of behaviour of the unreinforced mortar not (MN) and reinforced the 0.58% volume one of polypropylene fibres (MNP) in four-point bending test.

Figure 3 :

 3 Figure 3: Evolution of the four-point bending tensile strengths of the mortars with respect to the exposure temperature (rate: 2°C/min).

Figure 4 :

 4 Figure 4: Residual Young modulus (four-point bending test) vs. temperature of the heat treatment

Figure 5 :

 5 Figure 5: Evolution of the fracture energy with respect to the exposure temperature.

Figure 6 :

 6 Figure 6: Evolution of the stress intensity factor K I with respect to the exposure temperature.

Figure 7 :

 7 Figure 7 : Optical microscopy observations. (a) Calcined fibre. (b) Melted Fibre. (c) Partially melted fibre.

Figure 8 :

 8 Figure 8: a) mounting flame test, b) drilled external face,

Figure 9 :

 9 Figure 9: Comparison between ISO 834 fire and the experimental curve.

Figure 10 :

 10 Figure 10: Evolution of the temperature at different depths of the non-fibre mortar and fibre mortar.

Figure 11 :

 11 Figure 11: Thermal gradient within test specimens after one hour of thermal exposure during flame tests.

Figure 12 :

 12 Figure 12: Test bed to punching shear.

Figure 13 :

 13 Figure 13: Example of spalling of sample made with 0.58% polypropylene fibre and w/c=0.4.

Table 1 :

 1 Characteristics of test fibres (given SIKA)

						Young	Melting	Tensile
		Length	Diameter	Ratio	Density
	Fibre					Modulus	point	strength
		L (mm)	Ø (µm)	L/Ø	ρ (kg/m 3 )
						E (GPa)	(°C)	σ (GPa)
	Polypropylene	12	18	667	910	6	170	0,55

A Portland cement is used: CEM1 52,5 N CE CP2 NF. This cement is made up mainly of clinker 95%, whose details of the chemical and mineralogical compositions are reported in table 2.

Table 2 :

 2 Chemical and mineralogical composition of cement CEM I 52,5.

	Elements SiO 2 Al 2 O 3 Fe 2 O 3 CaO MgO K 2 O NaO 2 SO 3	IR	LOI Free CaO
	%	22.40 2.96	2.33	66.60 0.95 0.15 0.10	2.13 0.20 1.59 0.50
		C 3 S = 65.3	C 2 S = 18.6	C 3 A = 4.35	C 4 AF = 7.14
			IR : insoluble residue; LOI : loss on ignition.

Table 3 :

 3 Evolution of the bending tensile stress (MPa) with temperature.

	Mortar	T°C	20°C	400°C	500°C
	MN		7.3	3.4	2.0
	MNP		6.2	4.2	0.8

Table 4 :

 4 Young modulus (GPa) vs. temperature.

	Mortar	T°C	20°C	400°C	500°C
	MN		32.6	9.6	7.5
	MNP		31.6	13.8	3.8

Table 5 :

 5 Evolution of the fracture energy G f (J/m²) with respect to the exposure temperature.

	Mortar	T°C	20°C	400°C	500°C
	MN		1234	577	639
	MNP		3531	943	468

Table 6 :

 6 Evolution of the intensity factor K I(10 6 MPa.m ½ ) with respect to the exposure temperature.

	Mortar	T°C	20°C	400°C	500°C
	MN		6,3	2,4	2,2
	MNP		10,6	3,6	1,3

Table 7 :

 7 Residual punching shear effort of the plates after testing at flame.

		Non-heated	Heated
	Mortars	F max (kN)	F max (kN)
	Non-fibre	20.3	7.2
	Fibre	20.9	5.5


Which effective solar zenith angles to use for hourly irradiation measurements?

Philippe Blanc and Lucien Wald MINES ParisTech -PSL Research University, Sophia Antipolis cedex, France (lucien.wald@mines-paristech.fr) Time-series of in situ or satellite-derived measurements of hourly irradiation on horizontal surface are increasingly available. In many cases, such data are input to numerical procedures that can be performed only if a solar zenith angle can be associated to each measurement. Currently, only time stamp is given for each measurement; no angle is given. The solar zenith angle is greatly varying within 1 h and cannot be considered as constant. In this context, there is a practical request from companies, academics, researchers...: what is the best practice for computing this angle? The work presented here compares the performances of several existing practices and makes recommendations keeping in mind the practical aspects faced by practitioners, companies, academics, researchers... Highquality 1 min measurements of global, diffuse and direct collected by the BSRN station in Carpentras in France, are aggregated to yield hourly irradiation, and used to assess the performances. It is found that several practices may be used to compute in a satisfactory way the effective solar zenith angle. The most accurate results are obtained if one may obtain corresponding irradiation values that should be observed under clear-sky conditions. The McClear Service available at http://www.soda-pro.com is one Web service providing such estimates of cloud-free irradiation for several summarizations, from 1 min to 1 month. If one may not obtain such cloud-free estimates, then the most accurate results are obtained from using irradiation at the top of atmosphere. Performances show a tendency to decrease during sunrise and sunset hours.

  p e a n C o nfe re n c e


O Pitois
email: pitois@univ-mlv.fr

N Louvet

E Lorenceau

F Rouyer

Node contribution to the permeability of liquid foams

Keywords: Hydrodynamic resistance, drainage, permeability, Plateau border, vertex

This paper deals with the drainage of liquid foams. The liquid velocity is known to be related to viscous dissipation occuring within the elements of the liquid network, i.e. the channels and the nodes. When compared together, available values for the hydrodynamic resistance of a foam node appear to span over more than one order of magnitude. To clarify this point, we propose an alternative experimental method to estimate the value of this parameter. In contrast to previous experimental work performed on the foam scale, the node resistance is not treated as a fitting parameter, but instead it is measured directly on the microscopic scale. The results alow a consistent range of values to emerge for this parameter.

I Introduction

Recent experimental results [1] indicated that the permeability of a liquid foam results from both contributions of the liquid channels (Plateau borders) and their junctions (nodes). The permeability is dependent on the geometrical arrangement of the foam elements and also on the rheological properties of the interfaces that confine the liquid phase [1][2][3][4]. Two extreme regimes are expected as soon as the contribution of one element dominates the other: the channel dominated regime and the node dominated regime [5]. The transitions between these two regimes have been explained with the Boussinesq number -that compares the bulk dissipation with surface dissipation -as the only control parameter [1]. But, although the resistance of the channels has been found to fit this model for foam drainage, experimental [1][2][3][4] and numerical [6] studies do not allow a consistent value to emerge for the resistance of the nodes. More precisely, and as it is shown in the following, reported values span over more than one order of magnitude. The aim of this work is to clarify this point.

II Hydrodynamic resistance of foam nodes

The foam permeability is derived from the analysis of the flow on the microscopic scale, i.e. the foam network unit. For relatively dry foams, it is composed of 1 narrow Plateau border with radius of curvature r and ½ node. The total length of the unit cell L is related to the liquid volume fraction by: 𝜀 ≅ 𝛿 𝜀 𝑟 2 𝐿 2 ⁄ , with 0.171

   [2]
. Note that more complete expressions for  can be used to describe wetter foams, but as the dimensionless node resistance is expected to be independent on the liquid volume fraction, this latter relation is used in the following. Introducing the hydrodynamic resistances of the channel (R c ) and the node (R n ), the pressure difference over the network unit is related to the liquid flow rate Q with:

  2 cn P R R Q    (1)
The volume flow rate can be related to the average velocity through the foam channels v pb and the channel cross-sectional area

On the foam scale, the macroscopic pressure gradient G is related to the superficial velocity V S using Darcy's law:

S GV K   ( 3
)
where µ is the shear viscosity of the bulk liquid and K is the foam permeability (we recall that the superficial liquid velocity, i.e. the liquid flow rate flowing through the foam along a given direction divided by the cross-sectional area of the foam perpendicular to that direction, is 3

S pb

Vv   , where the factor 3 results from the angular average of the network units within the foam [2]). As we are solely interested in the resistance of the nodes, the resistance of the channels will be neglected in the following, i.e. 0 c R  . Thus, from ( 2) and ( 3), the foam permeability is given by:

2 3 2 12 2 3 an L K R     (4)
where

3 nn R r R  
is the dimensionless resistance. Several expressions are available for the foam permeability in literature. Comparisons of eq. ( 4) with these expressions provide several values for that can now be compared together (see table 1 for details). The values span over a large range although the foaming solutions used by the authors are similar. In fact, reported values in [2] and [1] are consistent, but it should be said that the corresponding authors used the same model for the foam permeability. In contrast, values reported in [6] and [4] are close together although the respective approaches are radically different. Considering the experimental data presented above, one can wonder if these large differences observed in the experimental values of the node resistance could be attributed to the assumptions made by the authors [1][2][3][4] about the network geometry. As the measured value for on the foam scale is model dependent -i.e. it is treated as a fitting parameter -the experimental study of a single node appears to be an interesting alternative method to estimate this value. The experimental setup used in this study has been previously described and more details can be found elsewhere [7,8]. We just recall that a Plateau border and the three adjoining films are formed on withdrawing a tripod from a reservoir containing the foaming solution, and that liquid can be injected through the channel at a controlled flow rate Q. For the present study, a gas bubble is released from the reservoir with a syringe allowing a node to be generated (see fig. 1). The total pressure difference P  can be measured over the system. P  results from the association in series of the upper Plateau border (pb1), the node (n), and the three lower Plateau borders (pb2), themselves in parallel association: (5) The determination of n P  requires the pressure drops 1 P  and 2 P  to be determined. These quantities are known to be dependent on the surface shear viscosity, µ S , characterizing the interface mobility. This latter parameter has been determined from the study of a single Plateau border, as reported in [7]. Then, the geometrical features of the channels (length and cross-sectional areas) are measured from image processing performed on pictures such as the one in figure 1. We considered that the junctions between the channels and the node correspond to the points of the profile where the curvature of the interface (within the vertical plane in figure 1) abruptly change or reverse. As the cross-sectional area slightly evolves along the PBs, the channels of length can be seen as composed of a set of N small borders, each with a length i l  , so that 0.0655 0.020 0.209

N ii i pbi i lB PQ AB              (6)
Note that the liquid flow rate through each channel pb2 is Q 2 = Q 1 /3. Besides, the effect of the longitudinal curvature of these channels is neglected in the estimation of the corresponding pressure drop.

For this study, TTAB (tetradecyltrimethylammoniumbromid 99 % purchased from Aldrich) was used in solution at a concentration C = 3 g/L (C/ CMC  3). This solution is similar to solutions used in [1][2][3][4]. With this solution, we found that µ S  3-4 10 -5 g/s and Bo  0.1, which is consistent with values reported in [1] for these parameters.

12 n

P P P P       

IV Results and discussion

In figure 2, the dimensionless parameter

   3 nn R r P Q  
is plotted as a function of the average liquid velocity through the upper Plateau border:

1 pb pb v Q A  (error bars correspond
to values of estimated with values of r at the inlet and the outlet of the node). For comparison, values reported in [1,2,3,4,6] are also presented: 𝑅 ̃𝑛 varies within the range 150-700 (see insert). This behaviour agrees qualitatively with the variation previously reported for a single Plateau border [7]. At low flow rates, pressure differences become close to the accuracy of the pressure measurement so that it is difficult to attribute the observed variation to an evolution of the mobility of the interface as this latter is sheared by the liquid flow. Note that this evolution can not be compared with other experimental data as the authors reported an average value of that has been fitted to their model over the range of investigated flow rates, and that such interfacial effects are not considered in numerical simulation [6]. At higher liquid velocities, full agreement is found with the values reported by Neethling et al. The average value of our data is very close to the numerical value proposed by Cox et al. This good agreement with [4] and [6] allows a consistent range of values to emerge: 150 <𝑅 ̃𝑛 < 400. In contrast, the full range of data reported in this paper does not overlap with the values reported in [1][2][3]. Could these discrepancies be explained by the influence of the Boussinesq number ? In fact, an estimation of this influence is available in [6]: in the limit of rigid interfaces, the contribution of the node is expressed as a length correction to the Plateau borders length. The corresponding value of is not explicitly reported, but from this paper, we could estimate 𝑅 ̃𝑛 = 4 × 0.03 × 49 𝛿 𝑎 2 ⁄ ≈ 225 . Comparing this value with the one corresponding to the mobile interface, i.e. 𝑅 ̃𝑛 ≈ 250 ∓ 25, it appears that although a slight increase is obtained for 𝑅 ̃𝑛 as a function of the interface mobility, this variation is not in quantitative agreement with figure 2, indicating that the observed discrepancies can not be explained by the effect of the interface mobility. As the values obtained for 𝑅 ̃𝑛 are model dependent in [1][2][3][4], a reasonable explanation for these discrepancies could be found in the network geometry assumed by the corresponding authors. A complete evaluation of this effect remains to be done. V Summary An alternative experimental method has been proposed to determine the hydrodynamic resistance of a foam node. In contrast to previous experimental evaluations, the node resistance is not treated as a fitting parameter, but instead it is directly measured on the microscopic scale. New results have been obtained and have been found to be in good agreement with one previous experimental work and with data obtained from the numerical study of this problem. This consistent set of values is not in accordance with the values recently proposed in the quantitative description of foam drainage, so that the present result seriously questions about the validity of this model, that overestimate the nodes resistance.

  At each element i is associated the cross-section pbi A , which is determined from the channel profile, and a local Boussinesq number

Figure 1 :

 1 Figure 1: Single foam node formed at the junction of four Plateau borders. The black vertical bar is one of the frame rods.

Figure 2 :

 2 Figure 2: Dimensionless resistance parameter of a foam node as a function of the average liquid velocity through the upper Plateau border (see figure 1). Error bars correspond to values of estimated with values of r at the inlet and the outlet of the node.

Table 1 :

 1 Comparison of available values for the node resistance 1 this expression has been deduced from[3] assuming that  < 1% and given by 𝜀 ≅ 𝛿 𝜀 𝑟 2 𝐿 2 ⁄

n R n R n R III Measurement of the resistance of a single node

Acknowledgements

We acknowledge the ANR financial support (ANR-05-JCJC-0234-01). n R


Henri Vallon

Andrea Cozza

Florian Monsef

Anne-Sophie Chauchat

Time-Reversed Excitation of Reveberation Chambers: Improving Efficiency and Reliability in the Generation of Radiated Stress

Keywords: Conversion efficiency, radiated test, reverberation chamber, time reversal, statistical electromagnetics

The ability of reverberation chambers to generate high-intensity field levels from relatively low-power input signals is reexamined for the case of time-reversed signals, proving that they lead to a higher efficiency. Moreover, the strong statistical spread typical of time-harmonic excitations can be dramatically reduced, thus improving the reliability of radiative tests, while limiting the need for a large number of independent realizations. The two excitation schemes are compared when forcing their respective input signals to display the same peak instantaneous power. Experimental results are provided, supporting the conclusions of our theoretical analysis.

I. INTRODUCTION

R EVERBERATION chambers (RCs) are capable of con- verting harmonic signals into high-intensity fields, by creating a pool of stored electromagnetic energy, thanks to their relatively long relaxation time constants. The stored energy can reasonably be expected to be uniformly spread over the entire volume of the RC, thus allowing for an ideally global excitation of an equipment under test (EUT).

Yet, this property also implies that, at each instant, only a portion of the energy injected into the RC interacts with the EUT. This observation is the starting point of our proposal for using time reversal (TR) as an alternative to harmonic or continuous-wave (CW) excitations. TR allows the generation of wavefronts focusing both in space and time; in other words, it allows concentrating a large portion of energy over a region of space at a given time. Studies of applications of TR to cavities are not new [START_REF] Draeger | One-channel time-reversal in chaotic cavities: Experimental results[END_REF]- [START_REF] Cozza | Statistics of the performance of time reversal in a lossy reverberating medium[END_REF], but they have been mainly concerned with demonstrating the possibility of producing focusing signals while propagating through non-anechoic media. We are rather interested in assessing how much TR can increase the peak intensity of the power used to stress an EUT.

This paper addresses this topic by analyzing the conversion efficiency of an RC when using TR signals rather than standard harmonic ones. The problem is tackled from a theoretical point of view, with the aim of providing simple predictive tools for the design of TR excitations. Our results prove that the benefits of TR come in three shapes: 1) ensuring a higher A. Cozza and F. Monsef are with PIEM, GeePs, UMR 8507, Centrale-Supelec, Univ Paris-Sud, UPMC, CNRS, 1 rue Joliot-Curie, 91192 Gif-sur-Yvette, France. A.-S. Chauchat is with THALES Communications & Security Gennevillers. H. Vallon is working with both institutions. Part of the work here reported was funded by the French National Research Agency through the grant ANR-12-ASTR-0005-01, MIMOCHIC project, within the framework of the 2012 ASTRID program. Contact e-mail: andrea.cozza@ieee.org field intensity, 2) reducing the uncertainty of the actual level of the generated field and 3) reducing the need for multiple realizations (i.e., stirring iterations) for the sake of producing a higher stress on the EUT. Even though not discussed in this paper, a fourth advantage should also be mentioned, namely the possibility of controlling the polarization of the resulting field, as proven in [START_REF] Cozza | Enforcing deterministic polarisation in a reverberating environment[END_REF]. The comparisons are carried out by acknowledging a fundamental limitation of power amplifiers, namely the maximum peak instantaneous power (PIP) they can handle. This constraint calls for an analysis differing from that presented in [START_REF] Cozza | Increasing peak-field generation efficiency of reverberation chamber[END_REF], while also taking into account the highly fluctuating behavior of RC-generated signals.

Only point-like focusing is considered in this paper, which means focal spots with spatial extensions subject to the diffraction limit; in practice, about half a wavelength wide. Extension to other kinds of wavefronts is possible by using generalized TR, a technique described in [START_REF]Emulating an anechoic environment in a wave-diffusive medium through an extended time-reversal approach[END_REF]. This same reference also discusses practical issues, e.g., how wavefronts are generated and moved around an EUT. This paper will only focus on how using TR excitations can improve the generation of electromagnetic stresses within an RC.

II. CONVERSION EFFICIENCIES

The point of view adopted in this paper is that an RC is a linear system that converts an input signal x(t), say a voltage, applied at the input port of an excitation antenna, into a vector electric field e(r, t) measured within the RC at the position r. In the rest of this paper we will drop the variable r, for the sake of brevity, and derive results that are valid for any r within the test volume, where we expect a diffuse-field approximation to hold. The average composite quality factor Q of the RC will be considered to be a good estimate of the composite quality factors experienced over the bandwidth used, around a working frequency f c . The average time constant of electromagnetic fields in the RC will thus be defined as τ = Q/πf c .

The field thus generated is expected to stress a portion of a circuit in an EUT, which will "receive" the power p r (t) as a linearly-polarized antenna, i.e.,

p r (t) = C r [e(t) • q] 2 = C r f 2 (t) (1)
where C r represents the sensitivity of the receiver, assumed to be non-dispersive, with no loss of generality, while q is the polarization-vector of the receiver. Similarly, we consider the input available power p i (t) = C i x 2 (t) applied at the antenna input port, with C i a constant of units equal to a conductance.

In view of the generation of a radiated stress, the conversion efficiency

η S = η[x S (t)] = p r (t) ∞ p i (t) ∞ = C r C i f 2 S (t) ∞ x 2 S (t) ∞ (2)
should be maximized, where the subscript S is used in order to distinguish different excitation signals. The rationale for choosing the L ∞ norm is twofold: the output PIP stressing the EUT should be as high as possible, while being generated by an input signal with the lowest PIP, in order to avoid any strain on power amplifiers and subsequent non-linear distortions or potential thermal stress. The q-aligned component of the electric field, f (t), is linearly related to the excitation signal through a transfer function Φ(ν), such that f (t) = F -1 {Φ(ν)X(ν)} (t), where F -1 {•} is the inverse Fourier transform between the frequency variable ν and the time variable t. The function X(ν) is the Fourier spectrum of the excitation signal. The following derivations will assume excitation signals with spectra centered at the frequency f c , covering the bandwidth

B T = f 2 -f 1 .
Harmonic signals are routinely used in standard operation of RC. In this case x CW (t) = A cos(2πf c t) and thus the conversion efficiency is

η CW = C r C i |Φ(f c )| 2 = C r C i W (f c ), (3)
where W (ν) = |Φ(ν)| 2 is a short-hand notation. In practice, the response of RCs is randomized by means of stirring techniques; as a result several realizations of Φ(ν) can be generated and the resulting efficiency can take very different values. This point will be discussed in sec. III.

In the case of time-reversed excitations, the input signal is chosen to be

x TR (t) = g(t) * φ(-t) = F -1 {G(ν)Φ ⋆ (ν)} (t), (4)
where G(ν) is the Fourier spectrum of the signal g(t) chosen to stress the receiver; * stands for the convolution operator applied to the time variable, while the superscript ⋆ represents the complex conjugate. Any kind of signal can be chosen for g(t), but time reversal has been demonstrated mainly with short pulses; a by-product of our derivation is that indeed better efficiencies are obtained by choosing a g(t) with the most compact support. Therefore, EUTs with a relaxation time longer than that of the cavity, typically in the microsecond range, are not suited to this technique, as their bandwidths would not cover enough degrees of freedom, as explained in the rest of the paper. This limitation can be overcome by using trains of short pulses; this kind of signals does not appear to have been considered so far in the literature and deserve investigations on their own as a practical solution.

Since the electric field f (t) at the receiver position will have a spectrum |Φ(ν)| 2 G(ν), if g(t) presents a maximum intensity at t = 0, the same will occur for f (t), given that |Φ(ν)| 2 is real-positive and corresponds to the spectrum of the auto-correlation function of φ(t), i.e., with its maximum at t = 0. Hence,

p r (t) ∞ = C r BT dν G(ν)W (ν) 2 = C r I 2 .
(

) 5
Solving the integral I by parts yields

I = G(f 2 )C W (f 2 ) - BT dνG ′ (ν)C W (ν), (6)
with G ′ (ν) the first derivative of G(ν) with respect to ν and

C W (ν) = ν f1 ds W (s) ≃ (ν -f 1 ) W (7)
a primitive of W (ν). The linear approximation in ( 7) assumes an average value W of W constant over B T ; indeed, typical RCs present slowly-varying statistical moments for the scalar components of the electric field, as shown in sec. IV. This approximation holds for ν -f 1 ≫ B c , where the B c = 1/τ is the coherence bandwidth of Φ(ν): in other words, as soon as the integration is carried out over enough significantly different portions of W (ν). Making use of this approximation in ( 6) and solving by parts the integral therein yields

I = W BT dν G(ν) = W ḠB T ( 8
)
and finally substituting it into (5) gives

p r (t) ∞ = C r W ḠB T . ( 9
)
This result holds for receivers modeled as in (1), i.e., presenting an flat frequency response over B T . An eventual constant group delay could be included, but is not explicited in (1) as it would not affect the PIP. When dealing with dispersive responses, ( 9) should be modified by weighting W (ν) by the transfer function of the receiver, before computing its average value over B T . Depending on its degree of frequency dispersion, this latter could reduce the frequency average and thus the peak received power. As thus, the case of a flat response must be regarded as a best case for TR excitations.

Before computing p i (t) ∞ , we need to recall that TR makes sense only for impulse-like signals g(t) covering a bandwidth B T ≫ B c , i.e., for B T /f c ≫ π/Q; under these conditions, g(t) has a time support T g ≪ τ , implying that the macroscopic features of φ(-t) ⋆ g(t) are dominated by those of φ(t).

Impulse responses in an RC can be modeled as random processes, with a root-mean-square (rms) envelope following

|φ(t)| 2 = A 2 o exp(-2t/τ ) [7, ch. 3] [8]
. The local deviations of the peak power of these kind of signals with respect to their rms envelope can be expressed as an overshoot factor

K = p i (t) ∞ /A o , ( 10
)
whose cumulative distribution function is

F K (x) = N k=0 erf xe k/BT τ / √ 2 , ( 11
)
with N ≥ 4B T τ ln 2, as demonstrated in [START_REF] Cozza | Stochastic Modeling of Impulse Responses Observed in Reverberating Environments[END_REF]. It is now sufficient to express A o as a function of

W (ν); since A 2 o = x 2 TR (0) , x 2 TR (0) = BT BT dν 1 dν 2 G(ν 1 )G ⋆ (ν 2 ) Φ(ν 1 )Φ ⋆ (ν 2 ) , ( 12
)
where

Φ(ν 1 )Φ ⋆ (ν 2 ) = |Φ(f c )| 2 µ(ν 2 -ν 1 ) (13)
is the frequency covariance function, which can be factorized into a slowly-varying power spectral density |Φ(ν)| 2 = W (ν) , here assumed to be constant over B T , such that W (ν) ≃ W , and a frequency-coherence function µ(ν 2 -ν 1 ) with a much narrower support, of the order of the modal bandwidth of the RC, as discussed in [START_REF]Multiple-Source Time-Reversal Transmissions in Random Media[END_REF]. For any signal g(t) such that its spectrum varies over scales much larger than B c , the coherence function operates as a sifting function

BT dν 2 G ⋆ (ν 2 )µ(ν 2 -ν 1 ) ≃ G ⋆ (ν 1 ) BT dν 2 µ(ν 2 -ν 1 ) = G ⋆ (ν 1 )B c , (14) yielding A 2 o = x 2 TR (0) = W B c BT dν |G(ν)| 2 , ( 15
)
where the last integral represents the energy of g(t), in the mathematical sense. Using ( 9) and ( 15), the conversion efficiency of time-reversed excitations is

η TR = C r C i B T B c W κ 2 K 2 , ( 16
)
with

κ 2 = B T Ḡ2 BT dν |G(ν)| 2 -1 ≤ 1 (17)
a shape factor attaining its maximum value for flat Fourier spectra.

III. POWER GAIN

How do η CW and η TR compare? In order to answer this question, two points should be considered: the random nature of W (ν) and the fact that, in case of harmonic excitations, an RC is typically operated jointly with a stirring technique. The assumption of an overmoded RC, as defined in [START_REF] Cozza | The role of losses in the definition of the overmoded condition for reverberation chambers and their statistics[END_REF], implies that W (ν) is bound to follow an exponential probability law [START_REF] Lehman | A statistical theory of electromagnetic fields in complex cavities[END_REF], [START_REF] Monsef | Variability and confidence intervals of the power measured in a reverberation chamber[END_REF]. As a result, η CW , when computed for a single realization, is widely dispersed around its mean value: e.g., η CW / η CW ∈ (0.05, 3.00) with 90 % probability. The poor predictability of the actual stress generated by an RC is improved by having N realizations, introducing stirring techniques [START_REF] Lehman | Characterization of the maximum test level in a reverberation chamber[END_REF]. This pragmatic approach is acknowledged by introducing

η N CW = max i∈[1,N ] η (i) CW (18)
i.e., the performance of the harmonic excitation will be evaluated as the highest one out of N realizations. As a result, x = η N CW / η CW will follow the probability law [START_REF] Ladbury | Evaluation of the NASA Langley Research Center mode-stirred chamber facility[END_REF], [START_REF] Hoijer | Maximum power available to stress onto the critical component in the equipment under test when performing a radiated susceptibility test in the reverberation chamber[END_REF] p

N (x) = N (1 -e -x ) N -1 e -x . ( 19
)
The strong skewness of this family of probability laws (see Fig. 1) implies that the statistical mode should be regarded as a more accurate estimator of typical values taken by η N CW , instead of its mean value. Therefore, the two conversion efficiencies will be compared by taking the ratio of their modes, a quantity hereafter referred to as the power gain

G N p = Mo [η TR ] Mo η N CW , ( 20
)
with Mo [•] the mode of a population of samples. The rationale for naming this quantity power gain is to be found in the possibility of obtaining the same order of magnitude of stress onto the EUT, while using a lower PIP at the input of the RC, i.e., at the output of the power amplifier. This observation is a direct consequence of having defined the conversion efficiency (2) as a quantity normalized to a fixed peak instantaneous input power.

As opposed to the case of harmonic excitation, time-reversal excitation is characterized by a self-averaging property [START_REF] Derode | Random multiple scattering of ultrasound. ii. Is time reversal a self-averaging process?[END_REF], as long as B T ≫ B c . This property is apparent in the definition of C W (ν) in [START_REF] Richardson | Reverberant Microwave Propagation[END_REF], where the integral is taken over a stationary random process W (ν) which takes only positive real values. The low coherence of W (ν) outside sub-bandwidths of the order of B c implies that the integral can be recast as a discrete sum

C W (ν) ≃ B c Nν i=1 W (s i ) (21) with s i -s i-1 = B c and ν -f 1 = N ν B c .
The stationarity of W (ν) over B T here translates into i.i.d. W (s i ); hence, the central-limit theorem allows stating that asymptotically C W , and ultimately f TR (t) ∞ , follow a normal distribution with relative standard deviation approximated by B c /B T . This property is confirmed in sec. IV. Use of a sufficiently large bandwidth can therefore be expected to reduce the uncertainty in the TR-generated field; this key point is discussed in secs. IV and V. Convergence to a weakly fluctuating f TR (t) ∞ invites the following approximation: to consider that the randomness of η TR is only due to x 2 TR (t) ∞ . Under this assumption, the pdf of η TR can be expressed as a function of the pdf of the random variable 1/K 2 , as apparent from ( 16) where p K (x) is the first derivative of F K (x) (to be computed numerically), and α = C r B T W κ 2 /C i B c . This pdf serves as a reference in sec. IV.

p K -2 (x) = (x/α) -3/2 2α p K α/K (22)
The mode, and also the mean of η TR are well approximated as

Mo [η TR ] ≃ η TR ≃ C r C i B T B c W (f c ) (κ/Mo [K]) 2 , ( 23
)
where

Mo [K] ≃ a ln b (4B T /B c ) (24)
was shown to be accurate to better than 1 % in [START_REF] Cozza | Stochastic Modeling of Impulse Responses Observed in Reverberating Environments[END_REF], when a = 0.749 and b = 0.678. For a harmonic excitation, the mode is

Mo η N CW = η CW ln N = C r C i W (f c ) ln N. ( 25
)
From ( 23) and (25), we can now express (20) as

G N p = B T B c κ 2 ln N a ln b (4B T /B c ) -2 . ( 26
)
It could seem surprising, at first sight, to see that for N = 1, G p is expected to be infinite. This is a direct consequence of the high probability of observing a weak field intensity for a single realization. In fact, the case N = 1 is of no practical use, as in any experimental setup typically involves at least a few tens of independent stirrer states.

Eq. ( 26) provides a direct estimate of how many realizations (e.g., stirrer steps) should be carried out in case of an harmonic excitation, in order to obtain a performance similar to that expected with TR excitations. Enforcing

G N p = 1, one would need N = exp B T B c a ln b (4B T /B c ) -2 (27)
independent realizations. A graphical representation of (27) in Fig. 2 shows that keeping up can quickly become unaffordable. Eq. ( 27) can also be interpreted as the ratio, between the bandwidth that should be spanned for CW frequency stirring, with respect to that needed by a TR signal in order to obtain the same efficiency. The underlying reason for this potentially huge difference in the performance is to be found in the collaborative nature of TR, where the N degrees of freedom are used coherently at the same time, as opposed to the CW case, where they are used separately in a non-coherent manner. Frequency (GHz) η CW (a.u.)

IV. EXPERIMENTAL VALIDATION

One of Supelec's RCs was used for the validation; its dimensions are 6.6 × 3.0 × 2.5 m 3 with an expected lowest usable frequency around 400 MHz. We did not study η N CW experimentally, since two decades of research on RCs have thoroughly validated the good accuracy of (19) when using harmonic excitations [START_REF] Lehman | Characterization of the maximum test level in a reverberation chamber[END_REF], [START_REF] Ladbury | Evaluation of the NASA Langley Research Center mode-stirred chamber facility[END_REF]. It is sufficient to estimate W (ν) , as discussed in sec. II. The validation therefore focused on the performance of TR excitation and the power gain, so that we did not apply any stirring technique, as TR is based on self-averaging. The presence of the stirrer was still useful for breaking symmetries in the cavity, as well as for reducing any direct illumination of the receiver from the source antenna.

A monocone antenna, mounted between a corner of the RC and the stirrer, was used for the excitation of the RC. Transfer functions between this antenna and a phase-sensitive electrooptical probe (Enprobe EFS-105) were measured by means of a vector network analyzer, collecting data over 5000 frequency samples in the bandwidth [1.8, 2.2] GHz. The probe was made to scan a hemispherical surface with a radius of 1 m, thanks to the low-perturbation robot depicted in Fig. 3, together with the probe. The residual field coherence, or correlation, as the robot moves was estimated to be about 80 %, thanks to its thin structure made out of glass fibers immersed in a plastic matrix. 3431 positions were thus considered, uniformly distributed with an average distance of 3.75 cm, i.e., λ/4 at 2 GHz. From this set of data we first computed η CW , shown in Fig. 4 to be indeed slowly varying. Field-related results were not corrected for the antenna factor of the probe, since our models refer to ratios of efficiencies affected by the same constants. Therefore efficiencies are all reported in arbitrary units (a.u.).

Time-domain responses φ(t) of the RC were computed at each position, for several bandwidths B T , up to 256 MHz. It was thus possible to estimate a time constant τ = 2.14 µs by means of a least-square fitting of φ 2 (t) ; τ is needed in order to compute the coherence bandwidth B c = 467 kHz required by our models.

Having direct access to the impulse responses, TR excitations can be defined as in (4), and their actual peak values recorded; all the results about TR excitations involve G(ν) = 1 over B T , i.e., a sine cardinale (sinc) pulse in the time domain. The histogram of the peak values, normalized to A o = W B c B T (see [START_REF] Hoijer | Maximum power available to stress onto the critical component in the equipment under test when performing a radiated susceptibility test in the reverberation chamber[END_REF]), is shown in Fig. 5, together with the theoretical pdf of the overshoot factor K, obtained by numerically differentiation of [START_REF] Lehman | A statistical theory of electromagnetic fields in complex cavities[END_REF]. These results refer to two bandwidths, B T = 32 and 256 MHz, corresponding to fractional bandwidths B T /f c equal to 1.6 % and 12.8 %, respectively. Both distributions present a good agreement, while their mode passes from 2.25 to 2.85, i.e., slightly more than a 25 % change for an eightfold increase in the excitation

bandwidth.

Conversely, the output field generated by these excitations is expected to be proportional to the bandwidth, as stated by [START_REF]Multiple-Source Time-Reversal Transmissions in Random Media[END_REF]. Fig. 6 supports this prediction for the mean value of TR output fields. In order to assess whether their mean value is a good representative of TR-generated fields, Fig. 7 presents the histograms obtained by computing ( f TR (t) ∞ -u)/ √ u, with u = (B T /B c ) W . If TR-generated fields were normally distributed according to the central-limit theorem, the previous operation should lead to a random variable well described by a standard normal distribution. The results in Fig. 7 confirms that this model is consistent for B T = 32 MHz, while for B T = 2 MHz the approximation, while not perfect, is still reasonable. This last observation is surprising, since for B T = 2 MHz, the input signal just covers slightly more than 4 coherence bandwidths, i.e., a very low number of degrees of freedom. These results have a strong effect on the confidence margin of TR-generated fields, as discussed in the next section.

The TR efficiency was then assessed, first for its mean value shown in Fig. 8. Its sub-linear dependence on B T is a direct consequence of how the maxima of RC impulse responses increase, although slowly, with B T [START_REF] Cozza | Stochastic Modeling of Impulse Responses Observed in Reverberating Environments[END_REF]. The slight discrepancy (about 6 %) between theory and experiments could seem surprising, since previous results were in very good agreement for both the numerator and the denominator of η TR . In fact, predicting the pdf of the ratio of random variables is not a simple matter, in particular in the case of interest, since the numerator and denominator are not independent, having a common origin in the impulse responses of the RC.

Having chosen to approximate the field generated by TR as a deterministic quantity is the reason for this residual error, which is confirmed in the comparison of the empirical and theoretical probability distributions (22) in Fig. 9. Finally, we present data confirming the validity of (26), shown in Fig. 10. Attention should be paid to the meaning of G N p ln N : it does not represent the final power gain, but a normalized version independent of N . The slight disagreement is a direct consequence of the one in the modeling of η TR .

V. CONFIDENCE INTERVALS

The good agreement of experimentally estimated TRgenerated peak fields with a standard Gaussian distribution implies that indeed output fields do converge as dictated by the central-limit theorem, i.e., they span a confidence interval C TR , e.g., with a 95 % probability, approximately equal to

C TR / η TR = 2 B c /B T , ( 28
)
to be compared with the one obtained for η

N CW C N CW / η N CW = ln 1 - N √ 0.025 1 - N √ 0.975 / ln N. ( 29
)
These confidence intervals, represented as functions of the available number of degrees of freedom, are shown in Fig. 11. The self-averaging property of TR excitations is the inherent cause for its faster convergence. It could therefore be disappointing to have η TR still featuring a pdf that does not appear to converge to a deterministic value. In fact, the reasons for the statistical dispersion of η TR are utterly different from those that are behind η CW . In the case of TR, the output field does converge to a deterministic value, as demonstrated in the previous section in Fig. 7; but when normalized to the input peak power, it is the stronger dispersion of this last quantity that increases the spread of η TR .

Conversely, in the case of CW excitations, the statistical dispersion only depends on the output fields, meaning that its actual value is hardly predictable. Therefore, in the case of TR excitations, a non deterministic η TR has a simple meaning: while the output field is predictable with a high degree of accuracy, the actual amount of peak input power is not. As such, it could require more, or less, peak power than expected on average, but this uncertainty will hardly have any major effect on the field intensity.

VI. CONCLUSIONS

The benefits of driving an RC with TR signals has been confirmed on three levels: 1) a higher conversion efficiency, 2) a reduced uncertainty in the actual field level, and this 3) without the need for multiple realizations (stirring). The increase in the efficiency may well be translated into a lightening in the peak-instantaneous power specifications for power amplifiers.

Our work has also proven that TR does not necessarily rhyme with wide-band signals, as acoustical TR has got us used to. Narrow bands can be envisaged as TR excitation signals, since the only constraint is B T /B c ≫ π/Q; with RCs easily displaying Q 1000, it should be possible to keep using standard power amplifiers while reaping TR benefits. Furthermore, a minimum frequency resolution can also be preserved.

The price to pay is the need for an overhaul in the use of RCs, as TR excitations require thinking in terms of spatiallyresolved stresses. A TR-driven RC works by focusing part of the injected energy over a reduced region of space during a short time interval. As a result, tests would need to be repeated when switching from one position to the next one, not unlike an anechoic chamber, where directions are tested one after the other. Depending on the size of the focal spot, and ultimately, on the desired spatial resolution and the EUT dimensions, the number of excitations required could vary wildly. Predicting this number is tricky, as a large number of potential implementations for TR-powered RC facilities can be devised. A discussion about the benefits and shortcomings of the various classes of TR implementations can be found in [START_REF]Emulating an anechoic environment in a wave-diffusive medium through an extended time-reversal approach[END_REF], where a special attention was given to the use of TRexcitations for test facilities.

Figure 1 .

 1 Figure 1. Probability distributions of η N CW / η CW for a varying number N of realizations.

Figure 2 .

 2 Figure 2. Equivalent number of independent realizations required in a CWdriven RC in order to ensure the same typical efficiency than a TR-driven one, as given by (27).

Figure 3 .

 3 Figure 3. The hemispherical positioner (a) and the electro-optical probe (b) used during the experimental validation.

Figure 4 .

 4 Figure 4. Estimate of η CW (grey dots) and its smoothed version (black line) obtained with a moving average over a 20 MHz bandwidth.

Figure 5 .Figure 6 .

 56 Figure 5. Empirical and theoretical pdfs of x TR (t) ∞ /Ao, for B T = 32 MHz (left) and 256 MHz (right).

Figure 7 .Figure 8 .

 78 Figure 7. Empirical and theoretical pdfs of ( f TR (t) ∞ -u)/ √ u, with u = (B T /Bc) W , for B T = 2 MHz (left) and 32 MHz (right).

Figure 9 .Figure 10 .

 910 Figure 9. Empirical pdfs of η TR /(B T W /Bc) and approximate pdf obtained by assuming a deterministic output, as defined in (22) for: B T = 32 MHz (left) and B T = 256 MHz (right).

Figure 11 .

 11 Figure 11. Relative confidence intervals of the conversion efficiencies for a 95 % probability, as a function of the number of degrees of freedom.


Matthias Caenepeel
email: mcaenepe@vub.ac.be

Fabien Seyfert
email: fabien.seyfert@inria.fr

Yves Rolain

Martine Olivi

Parametric Modeling of the Coupling Parameters of Planar Coupled-Resonator Microwave Filters

The design of planar coupled-resonator microwave filters is widely based on coupling matrix theory. In this framework a coupling matrix is first obtained during the synthesis step. Next this coupling matrix is physically implemented by correctly dimensioning the geometrical parameters of the filter. The implementation step is carried out using simplified empirical design curves relating the coupling coefficients and geometrical parameters. The curves typically only provide initial values and EM optimization is often needed such that the filter response meets the specifications. This paper proposes to extract parametric models that relate the filters design parameters directly to the coupling parameters. The advantage of such models is that they allow to tune the filter in a numerically cheap way and that they provide physical insight in the filters behavior. This paper explains how such models can be extracted from EM simulations and illustrates the technique for the design of an 8 pole cascaded quadruplet filter in a microstrip technology.

I. INTRODUCTION

Coupling matrix theory is widely used to design narrowband bandpass microwave filters [START_REF] Cameron | General coupling matrix synthesis methods for chebyshev filtering functions[END_REF], [START_REF] Cameron | Microwave filters for communication systems[END_REF] in which the filter is modeled as a low-pass coupled resonator circuit (Fig. 1) [START_REF] Atia | New types of waveguide bandpass filters for satellite transponders[END_REF]. The design process begins therefore by the derivation of a filtering characteristic meeting the filtering specifications. In a second step this characteristic is realized by a circuit. Next the circuit elements (coupling parameters) must be implemented physically by correctly dimensioning the geometrical parameters of the filter. An initial dimensioning of the filter is carried out by means of simplified empirical design curves relating the coupling parameters and geometrical parameters of the filter [START_REF] Hong | Microstrip filters for RF/microwave applications[END_REF]. These curves however typically do not take into account more complex interactions, such as loading effects of the resonators by the couplings, which require that the initial values must be optimized to meet the filter specifications. Several electomagnetic (EM) optimization techniques exist in the literature [START_REF] Cameron | Microwave filters for communication systems[END_REF]. Although these techniques prove to be successful, there are still some disadvantages. The computation time grows rapidly with the complexity of the design (number of geometrical parameters). When the filter specifications change even slightly, the whole optimization process must be relaunched. Moreover the optimization process does not provide any physical insight in the filters behavior. Literature shows that the coupling parameters are smooth functions of the controlling geometrical parameters [START_REF] Hong | Microstrip filters for RF/microwave applications[END_REF], [START_REF] Amari | Space-mapping optimization of planar coupled-resonator microwave filters[END_REF]. Therefore we propose to approximate this relation as a quadratic multivariate polynomial. The multivariate character stems from the fact that we also take into account several parameters to model second order effects such as the loading of the resonators. The main advantage is that these models are numerically cheap to evaluate and provide physical insight in the filters behavior. Moreover they can be re-used in various design scenarios. This paper explains how such parametric models can be extracted from EM simulations. A crucial step in the modeling process, is the extraction of the coupling parameters from the simulated scattering data (Sparameters).

State-of-the-art gradient-based parameter extraction methods optimize the generalized low-pass network such that the networks frequency response meets the measured or simulated response [START_REF] Harscher | Automated filter tuning using generalized low-pass prototype networks and gradient-based parameter extraction[END_REF], [START_REF] Harscher | Em-simulator based parameter extraction and optimization technique for microwave and millimeter wave filters[END_REF]. In the case where multiple solutions are possible, these methods do not necessarily converge to the implemented circuit. Section II shortly describes a three stage extraction process that overcomes this drawback. The Sparameters and network parameters are related as follows [START_REF] Olivi | Identification of microwave filters by analytic and rational h2 approximation[END_REF]

S(s) = I + C(sI -A) -1 C t with C = i √ 2R in 0 . . . 0 0 . . . 0 i √ 2R out , (1)
A = -R -iM. (2)
The first step identifies a rational matrix from the simulated S-parameters [START_REF] Olivi | Identification of microwave filters by analytic and rational h2 approximation[END_REF]. The second step finds all possible circuital realizations of the previously computed rational approximation with a specified coupling topology. The third step finally deciphers the physically implemented matrix. Section III describes how, given the physical coupling parameters, the quadratic multivariate polynomials are estimated in a least-squares sense. Eventually the modeling technique is illustrated for a design of an 8 pole cascaded quadruplet filter in microstrip technology (Fig. 2). The design example shows that the models can be used to improve the filter response drastically.

II. COUPLING PARAMETER EXTRACTION

This section describes the three stage parameter extraction procedure. First a canonical solution is extracted. Next all possible solutions corresponding to the desired topology are identified. Finally the physically implemented coupling matrix is identified among these solutions using prior knowledge. A more detailed explanation can be found in [START_REF] Caenepeel | Microwave filter design based on coupling topologies with multiple solutions[END_REF]. The first step identifies a rational matrix from the simulated Fig. 1: Low-pass prototype network model for a general crosscoupled resonator filter.

S-parameters.

The scattering data is first transformed to the low-pass domain. Next a rational common denominator matrix is identified using the toolbox PRESTO-HF [START_REF]PRESTO-HF : a matlab toolbox dedicated to the identification problem of low pass coupling parameters of band pass microwave filters[END_REF] . The order of the common denominator n is chosen to be equal to the number of resonators. The rational matrix is then transformed into the canonical arrow form of the coupling matrix. The second step finds all possible realizations for the given rational matrix corresponding to a given topology. This topology corresponds to the topology used during the synthesis step. It is however possible that for certain geometrical values, the corresponding S-parameters can not be realized with the ideal topology. E.g. the response is asymmetric while the topology is only able to realize symmetric responses. Such a situation might occur when frequency offsets and/or extra couplings are present in the physical structure. To handle these situations the allowable topology is extended to a topology that is close to the original topology but still allows to capture these effects. This corresponds to expanding the class of realizable responses of the filter. Next the arrow form obtained during step one is transformed to the new topology using the toolbox DEDALE-HF [START_REF]DEDALE-HF : a matlab toolbox dedicated to the equivalent network synthesis for microwave filters[END_REF]. The only drawback of this extended topology is that it concentrates the extra couplings at fixed positions in the matrix, while physically they might occur somewhere else. To compensate for this unwanted effect, the solutions found by DEDALE-HF are optimized by means of similarity transformations on the coupling matrix. This optimization process minimizes the influence of the extra couplings by redistributing them over the whole matrix. The influence of the extra couplings is expressed by a relative measure c:

c = Σ(M extra ) 2 Σ(M all ) 2 (3)
which is the sum of squares of the extra couplings (M extra ) over the sum of the squares of all the couplings (M all ). During the last step of the process, the physically implemented coupling matrix is estimated among the solutions found during step 2. The choice of this matrix is based on some physical assumptions. First of all we assume that the couplings present in the original topology are dominant in the implemented filter. Therefore only those solutions for which c is sufficiently small are considered. Since we use the extraction method in a design context, we assume that the implemented couplings are close to the ideal ones found during the synthesis step. Therefore the solutions for which the 2-norm of the difference to the ideal coupling matrix is minimal is chosen. To model the couplings as a function of the geometrical parameters, the filter is simulated for several geometrical values. This information allows to check whether the choice of the physical matrix is consistent with the geometry. E.g. when the spacing between resonator i and k is varied, this should mainly affect the coupling M ik . If this is not the case, the solution that is coherent with the geometric variation, is chosen instead.

III. MULTIVARIATE QUADRATIC APPROXIMATION

The design curves that are available in the literature show that the behavior of coupling parameters as a function of the geometrical parameters is smooth [START_REF] Hong | Microstrip filters for RF/microwave applications[END_REF], [START_REF] Amari | Space-mapping optimization of planar coupled-resonator microwave filters[END_REF]. Therefore we propose to model the coupling parameters using multivariate polynomials of a maximal degree of 2. Since the process described in section II, extracts the physical coupling matrix, it is possible to relate a coupling parameter directly to a set of geometrical parameters. We illustrate this with the example of 2 coupled resonators. It is known that the coupling between 2 resonators depends strongly on the spacing between the resonators, but is also affected by a possible offset introduced by spacings between other adjacent resonators (Fig. 3). In this case we propose to use the following model:

M ik = a 1 s 2 ik + b 1 s ik + b 2 d + c 1 (4)
where s ik is the spacing between the resonators and d represents the offset. In this case 4 coefficients a 1 , b 1 , b 2 and c 1 must be estimated. This requires at least 4 different EM simulations. The mesh that is used for the different simulations however might differ slightly since it depends on the geometry. This difference together with the error norm between the rational model and the EM simulation used during the extraction step, introduce an error on the extracted coupling matrices. To avoid that the parametric model of the coupling parameter would model these effects as well, we propose to use at least 2 times more simulations than what is minimally needed to estimate the coefficients. The geometrical values for which the filter must be simulated are selected as follows. First the initial values for the geometrical parameters are determined using the design curves. Next the physical coupling matrix for this structure is extracted. This coupling matrix gives an idea of the distance between the initial design and the ideal design. This allows to determine a validity interval around each design parameter. Values for each geometrical parameter are then randomly generated within this interval. The number of values corresponds to the proposed number of simulations needed to estimate the parametric models. The filter is simulated for all these sets of values. Once the EM simulations are gathered, the coupling parameters are extracted using the method described in section II. The couplings together with the geometrical parameters are then used to estimate a parametric model for each coupling parameter in least-squares sense.

IV. EXAMPLE

The design of an 8-pole cascaded quadruplet microstrip filter consisting of square open-loop resonators (Fig. 2) is used to illustrate the proposed method. The filter is designed for a center frequency f c of 1 GHz and fractional bandwidth F BW of 0.06. The ideal coupling matrix is synthesized using DEDALE-HF. The synthesis step yields 2 possible solutions among which one is chosen. The filter is implemented in a RT/duroid substrate with a relative permittivity 10.2. The initial dimensions of the filter are obtained using empirical formulae derived for the coupling coefficients between square open-loop resonators [START_REF] Hong | Couplings of microstrip square open-loop resonators for cross-coupled planar microwave filters[END_REF]. The filter is simulated using ADS Momentum 2014 [START_REF]Advanced design system[END_REF]. Fig. 4 and5 show that the responses can clearly be improved with respect to center frequency, insertion loss and bandwidth.

There are three types of coupling parameters that potentially need to be tuned and thus modeled: • Inter-resonator couplings M ik where i = k

• Center frequency offsets M kk • Normalized input/output impedances R 1 and R 2 re- spectively
It is well known that inter-resonator couplings mainly depend on the spacing between the resonators. Spacings between other resonators can introduce an extra offset (Fig. 3). E.g. the coupling between resonator 1 and 2 mainly depends on the distance s 12 . The spacing between resonator 1 and 4, labeled s 14 and the spacing between 2 and 3, labeled s 23 introduce an offset |s14-s23|

2

. Thus we can model this coupling as:

M 12 = a 1 s 2 12 + b 1 s 12 + b 2 ( |s 14 -s 23 | 2 ) + c 1 (5)
The same reasoning can be repeated for the other interresonator couplings.

The diagonal elements M kk express the difference between the center frequency of the filter and the resonance frequency f k of resonator k. f k is directly related to the length of the resonator and thus to g k , but also depends on the loading of the resonator influences f k . This loading is modeled by means of the distances between the neighboring resonators. Therefore we propose as a model:

M kk = a 1 g 2 k + b 1 g k + b 2 s k-1,k + b 3 s k,k+1 + b 1 (6)
if k = 1, n and

M kk = a 1 g 2 k +b 1 g k +b 2 s k-1,k +b 3 s k,k+1 +b 4 t in/out +c 1 (7)
if k = 1 or n. The normalized input and output impedances mainly depend on the position of the input and output feeding lines. The length of the lines however also has an effect. Therefore we propose a model:

1 = a 1 t 2 in + b 1 t in + b 2 g 1 + c 1 (8) R 2 = a 1 t 2 out + b 1 t out + b 2 g n + c 1 (9)
This reasoning shows that the most complicated model requires the estimation of 6 coefficients. We propose to perform 14 EM simulations to estimate them. Since the coupling parameters are not too far from the ideal ones, we propose to simulate the filters for values in an intervals varying from 0.3 up to 0.5 mm for the design parameters s ik , since the offset from f c are also not too large we use the same intervals for g k . For t in and t out we take an interval of 0.5 mm.

After the EM simulations are gathered, the coupling matrices are extracted using the method described in section II. Fig. 6 shows the extracted inter-resonator couplings for each simulation. Eventually the simulations are used to extract the parametric models given in equations 5-9. Fig. 7 shows the result of the model for M 12 for all of the simulations. It shows that there is model error, which is expected due to simulation meshing and parameter extraction errors.

Once the parametric models are extracted they can be used to fine-tune the filter. The models form a set of non-linear equations. This set is solved in Matlab using the function fminunc. The filter is next re-simulated for the found geometrical parameters. The result is shown in Fig. 4 and5. The response can still be improved, but remark that only 14 EM simulations were required to tune a filter with 19 design parameters.

V. CONCLUSION

This paper introduces parametric models that relate the coupling matrix parameters to the design parameters of the filter. The advantage of the models is that they are cheap to evaluate and allow to improve the filters response in a numerically cheap way. We explain how the models can be extracted using EM simulations. Moreover we illustrate the utility of the models with a design example. The design example shows that the models allow to improve the filters response drastically using only 14 EM simulations.

Fig. 2 :

 2 Fig. 2: Top-view of a square open-loop resonator cascaded quadruplet filter and the design parameters s ik , g i , t in and t out .

Fig. 3 :Fig. 4 :

 34 Fig.3: Top-view of 2 coupled resonators, where s ik is the spacing between them and d is a possible offset.

Fig. 5 :

 5 Fig. 5: Magnitude of S 11 for the initial (red) and final (blue) design.

Fig. 6 :Fig. 7 :

 67 Fig. 6: The extracted inter-resonator couplings for each simulation.

R

This work is sponsored by the Vrije Universiteit Brussel, dept. ELEC, Pleinlaan 2, 1050 Brussel Belgium, INRIA Sophia-Antipolis, Fund for Scientific Research (FWO-Vlaanderen), the Flemish Government (Methusalem), the Belgian Government (IUAP-VII) and the Strategic Research Program of the VUB (SRP-19).