Léa Penazzi
email: lea.penazzi@univ-lorraine.fr

Yves Jannot

Johann Meulemans

Olivier Farges

Vincent Schick

Influence of radiation heat transfer on parallel hot-wire thermal conductivity measurements of semi-transparent materials at high temperature

Keywords:

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

As previously pointed out in [START_REF] Jannot | An improved model for the parallel hot wire: Application to thermal conductivity measurement of low density insulating materials at high temperature[END_REF], knowledge of the thermal conductivity of high temperature insulating materials is of great importance for the control of industrial processes. It has also been stated that the parallel hot-wire method is one of the most suitable method to measure the thermal conductivity at high temperature [START_REF] Jannot | An improved model for the parallel hot wire: Application to thermal conductivity measurement of low density insulating materials at high temperature[END_REF][START_REF] Jannot | Apparent thermal conductivity measurement of anisotropic insulating materials at high temperature by the parallel hot-wire method[END_REF], moreover it is a standardized method [START_REF]Refractory Materials-Determination of Thermal Conductivity-Part 2: Hot-Wire Method (Parallel), Standard[END_REF]. Some recent developments in the modeling of this method have led to a better accuracy of the measurement [START_REF] Jannot | An improved model for the parallel hot wire: Application to thermal conductivity measurement of low density insulating materials at high temperature[END_REF][START_REF]Refractory Materials-Determination of Thermal Conductivity-Part 2: Hot-Wire Method (Parallel), Standard[END_REF]. Nevertheless, the estimation process is more often based on a purely conductive model, considering an equivalent thermal conductivity that is the sum of the thermal conductivity due to conduction and of the radiative thermal conductivity [START_REF] Jannot | An improved model for the parallel hot wire: Application to thermal conductivity measurement of low density insulating materials at high temperature[END_REF][START_REF] Jannot | Apparent thermal conductivity measurement of anisotropic insulating materials at high temperature by the parallel hot-wire method[END_REF][START_REF]Refractory Materials-Determination of Thermal Conductivity-Part 2: Hot-Wire Method (Parallel), Standard[END_REF][START_REF] Zhang | Simultaneous measurement of thermal conductivity and thermal diffusivity of solids by the parallel-wire method[END_REF][START_REF] Grazzini | Measuring thermal properties with the parallel wire method: A comparison of mathematical models[END_REF][START_REF] Santos | Numerical and experimental determination of minimum and maximum measuring times for the hot wire parallel technique[END_REF]. This approximation holds true for optically thick media.

Some authors investigated the effect of radiation heat transfer without making this assumption. Ebert and Fricke [START_REF] Ebert | Influence of radiative transport on hot-wire thermal conductivity measurements[END_REF] simulated the heat transfer in a hot-wire measurement with a one dimensional coupled conduction-radiation model valid only for optically thick media. The radiation emitted by the heating wire was also neglected. Their theoretical study showed that the hot-wire method is applicable only if the extinction coefficient exceeds a critical value. Gross and Tran [START_REF] Gross | Radiation effects on transient hot-wire measurements in absorbing and emitting porous media[END_REF] simulated the heat transfer in a hot-wire measurement with a one dimensional coupled conduction-radiation model solved by the discrete ordinates method. They showed that this critical value increases with the temperature. The influence of the proportion of diffusion/absorption in the extinction coefficient was not addressed in their theoretical study. Coquard et al. [START_REF] Coquard | Experimental and theoretical study of the hot-wire method applied to low-density thermal insulators[END_REF] developed a two dimensional coupled conduction-radiation model solved by the discrete ordinates method. They used it to study the accuracy of the estimation of the thermal conductivity by the hot-wire method, they highlighted the estimation bias for materials of low density or low optical thickness. Nevertheless, the study was applied to only one material (an EPS foam) and at ambient temperature. Daouas et al. [START_REF] Daouas | Solution of a coupled inverse heat conduction-radiation problem for the study of radiation effects on the transient hot wire measurements[END_REF] simulated the heat transfer in a parallel hot wire measurement with a one dimensional coupled conduction-radiation model solved by the discrete ordinates method. A sensitivity study showed the possibility of estimating the absorption coefficient in addition to the thermal parameters. However, the study was performed at room temperature on a very weakly https://doi.org/10.1016/j.ijthermalsci. absorbing and non-scattering material. Zhang et al. [START_REF] Zhang | Numerical study of the influence of thermal radiation on measuring semi-transparent thermal insulation material with hot wire method[END_REF] performed a numerical study simulating the coupled conduction-radiation heat transfer in a purely absorbing medium in both cases of hot-wire and parallel hot-wire measurements with discrete ordinates method. The numerical analysis shows that the thermal conductivity of insulated material with low extinction coefficient measured by hot-wire device is overestimated at high temperature. However, the study was conducted on only one material with unique density and thermal conductivity. The chosen thermal conductivity value was very low (𝜆 = 0.02 W m -1 K -1 ) and the maximum temperature does not exceed 1000 °Celsius.

The purpose of this study is to establish the conditions under which the temperature evolution in a parallel hot-wire measurement can be correctly represented by a pure conductive model using an effective thermal conductivity. A theoretical study using the P1 approximation method implemented in COMSOL Multiphysics ® will be performed for temperatures ranging from ambient to 1500 °C, for both absorbing and scattering media and for several densities, extending the range of materials and temperatures of the previous studies. When the purely conductive model applies, a new estimation method is proposed that improves both thermal conductivity and thermal capacity estimations. A calibration process of the distance between the hot wire and the thermocouple will also be proposed, enabling a more accurate estimation of the volumetric heat capacity. The results will be supported by an experimental study carried out on three insulating materials with densities ranging from 581 kg m -3 to 910 kg m -3 and at temperatures ranging from 20 °C to 1200 °C. Finally, a method enabling the estimation of the extinction coefficient from thermal conductivity measurements at various temperatures is presented and successfully applied to the three tested materials.

Principle and models

Principle

The schematic diagram of the classical parallel hot-wire set-up is presented in Fig. 1. A resistive heating wire with a radius 𝑟 w is inserted in a groove placed on the surface of the bottom sample. A thermocouple is inserted in another groove placed at the surface of the same sample at distance 𝑑 of the heating wire. A second sample with the same dimensions as the bottom sample is then placed on top of it. A constant heat flow rate is then imposed in the heating wire by passing an electrical current with a constant intensity 𝐼. The equations of the system are:

⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ 𝜌𝑐 𝜕𝑇 𝜕𝑟𝑡 = 𝛁 ⋅ (𝜆 c 𝛁𝑇 ) (a) 𝑇 (𝑟, 0) = 0 (b) lim 𝑟→∞ 𝑇 (𝑟, 0) = 0 (c) 𝜑 𝑤 = 𝜌 w 𝑐 w 𝜋𝑟 2 w 𝜕𝑇 w 𝜕𝑡 + 𝑇 w -𝑇 (𝑟 w ) 𝑅𝑙 w (d) (1)
A complete model of the system can be found in [START_REF] Jannot | An improved model for the parallel hot wire: Application to thermal conductivity measurement of low density insulating materials at high temperature[END_REF], it is based on the following hypotheses:

• The samples are infinite in the 𝑂𝑥 and 𝑂𝑧 directions (semi-infinite medium hypothesis).

• There is no temperature gradient along the 𝑂𝑦 direction (''long wire'').

• The heat flow rate 𝜑 w produced in the heating wire is constant.

• The mass of the thermocouple is negligible.

• Radiation heat transfer is negligible or may be represented by an effective thermal conductivity.

Jannot and Degiovanni [START_REF] Jannot | An improved model for the parallel hot wire: Application to thermal conductivity measurement of low density insulating materials at high temperature[END_REF] established the following expression for the Laplace transform of the thermocouple temperature:

𝜃(𝑑, 𝑝) = 𝜑 w 𝑝𝑙 w 𝐾 0 ( √ 𝑝 𝑎 𝑑) 𝜌 w 𝑐 w 𝜋𝑟 2 w 𝑝𝐾 0 (√ 𝑝 𝑎 𝑟 w ) + 2𝜋𝜆 c √ 𝑝 𝑎 𝑟 w [ 1 + 𝜌 w 𝑐 w 𝜋𝑟 w 2 𝑝𝑅𝑙 w ] 𝐾 1 (√ 𝑝 𝑎 𝑟 w ) (2)
The temperature rise 𝑇 (𝑑, 𝑡) is calculated by applying an inverse Laplace transform to (2) using the De Hoog algorithm [START_REF] De Hoog | An improved method for numerical inversion of Laplace transforms[END_REF].

General equations

It is assumed here that the medium is gray. The coupled conduction and radiation heat transfer in transient hot-wire measurements is described by the following energy balance:

𝜌𝑐 𝜕𝑇 𝜕𝑡 = -𝛁 ⋅ (𝜙 𝜙 𝜙 c + 𝜙 𝜙 𝜙 r ) (3)
Where 𝜙 𝜙 𝜙 c is the conductive heat flux given by ( 4) and 𝜙 𝜙 𝜙 r is the radiative flux given by [START_REF] Zhang | Simultaneous measurement of thermal conductivity and thermal diffusivity of solids by the parallel-wire method[END_REF].

{ 𝜙 c = -𝜆 c 𝜕𝑇 𝜕r (a) 𝛁 ⋅ 𝜙 𝜙 𝜙 r = 4𝑘 a 𝑛 2 𝜎𝑇 4 -𝑘 a ∫ 4𝜋  (𝝎) d𝛺 (b) (4)
(𝝎) is the radiant intensity in the 𝝎 direction given by the radiative transfer equation (RTE) (Eq. ( 5)):

𝜕 (𝝎) 𝜕𝑠 = - ( 𝑘 a + 𝑘 s )  (𝝎) + 𝑘 a 𝑛 2 𝜎𝑇 4 + 𝑘 s 4𝜋 ∫ 4𝜋  ( 𝝎 ′ ) 𝑃 ( 𝝎 ′ → 𝝎 ) d𝛺 ′ (5)
Where 𝑘 a is the absorption coefficient, 𝑘 s , the scattering coefficient, 𝑛, the refractive index and 𝜎 is the Stefan-Boltzmann constant.

Model M1: P1 approximation

The P1 approximation, corresponding to a spherical harmonics method (PN) truncated at the first order, enables solving the RTE presented in [START_REF] Grazzini | Measuring thermal properties with the parallel wire method: A comparison of mathematical models[END_REF]. The P1 approximation relies on the assumption that the scattering is linear isotropic [START_REF] Michael | Chapter 16 -the method of spherical harmonics (PNapproximation)[END_REF]. This assumption is justified in this present configuration as the medium is optically thick, 𝜏 ≫ 1, where 𝜏 is the optical thickness defined by the integral of the extinction coefficient 𝛽 along a typical optical path 𝜏 = ∫ 𝑠 0 𝛽d𝑠. From a computational point of view this approximation has a limited impact because it introduces only one additional degree of freedom for the incident radiation 𝐺 (W m -2 ), which is a scalar quantity and adds a heat source or sink to the temperature equation to account for radiative heat transfer contributions. This method, however, fails to accurately represent cases where the radiative intensity propagation dominates over its diffusivity or where the scattering effects cannot be described by a linear isotropic phase function. The P1 approximation is a simplified approach, which is expected to be fairly good for absorbing and highly scattering media at large optical distances from boundaries or interfaces that have a strong variation of temperature and radiative characteristics of the medium [START_REF] Howell | Thermal Radiation Heat Transfer[END_REF]. The P1 approximation solves the RTE equation, expressed in [START_REF] Grazzini | Measuring thermal properties with the parallel wire method: A comparison of mathematical models[END_REF], by solving the following equation:

𝛁 ⋅ (-𝐷 𝑃 1 𝛁𝐺) = -𝑄 r ( 6
)
where 𝐺 is the incident radiation (Eq. ( 7)), 𝐷 𝑃 1 is the P1 diffusion coefficient in the case of an isotropic diffusion (Eq. ( 7)) and 𝑄 r is the radiative volumetric heat source (Eq. ( 7)).

⎧ ⎪ ⎨ ⎪ ⎩ 𝐺 = ∫ 4𝜋 0  (𝝎) d𝛺 (a) 𝐷 𝑃 1 = 1 3𝑘 a + 3𝑘 s (b) 𝑄 r = 𝛁 ⋅ 𝜙 𝜙 𝜙 r = 𝑘 a (𝐺 -4𝑘 a 𝑛 2 𝜎𝑇 4 ) (c) (7)
Regarding the boundary conditions, at the hot wire surface, thermocouples surfaces and sample external surfaces are considered as opaque surfaces. The following equation accounts for the radiation emitted/absorbed by the boundaries and is known as the Marshak's boundary condition [START_REF] Marshak | Note on the spherical harmonic method as applied to the milne problem for a sphere[END_REF]:

𝒏 ⋅ ( -𝐷 𝑃 1 𝛁𝐺 ) = 𝜙 r,net (8)
Where 𝜙 r,net is the net radiative heat flux at the boundary and 𝒏 the outward unit direction on the surface. In the following simulations, surfaces of the hot wire, thermocouples and the external surface of the sample were considered gray with an approximate common emissivity value of 𝜀 = 0.8. The net radiative heat flux is then expressed as follows:

𝜙 r,net = 𝜀 2(2 -𝜀) (4𝑛 2 𝜎𝑇 4 -𝐺) (9)

Model M2: Rosseland approximation

If the medium is gray and optically thick, the radiative flux divergence is expressed as follows:

𝜙 r = -𝜆 r 𝜕𝑇 𝜕𝑟 ( 10
)
With the expression of the radiative thermal conductivity [START_REF] Jannot | Thermal Properties Measurement of Materials[END_REF]:

𝜆 r (𝑇 ) = 16 3 𝑛 2 𝜎 𝑇 3 𝛽 ( 11
)
where 𝛽 is the extinction coefficient of the medium given by 𝛽 = 𝑘 a + 𝑘 s . We are thus brought back to a pure conduction problem with an effective conductivity 𝜆 = 𝜆 c + 𝜆 r .

Materials and experimental device

An experimental study was carried out to confirm the theoretical results. Measurements were undertaken on an isotropic material: SilPower ® (Saint-Gobain Quartz) with three different densities: 581, 741 and 910 kg m -3 referenced respectively as SG, MG and HG. SilPower ® Rigid Silica is a high performance infra-red reflector obtained from Quartzel ® fibers having an SiO 2 content ≥ 99.95 % (Saint-Gobain, 2021).

The dimensions of the samples are 230 × 120 × 60 mm 3 as recommended by the standard ISO 8894-2:2007. The parallel hot-wire (PHW) measurements were conducted at 20, 200, 400, 600, 800, 1000 and 1200°C. Three experiments were carried out for each temperature and for each density. Fig. 2 presents a schematic view of the measurement set-up. The heating wire made of Nickel-Chrome 80/20 has a diameter of 0.5 mm and it is inserted in a groove with a 0.5 × 0.5mm 2 section. The groove with the same length as the sample was machined in the middle of one face of the sample.

To study the influence of the distance between the heating wire and the thermocouple, the temperature rise was measured at two locations, respectively at 𝑑 1 ≈10 mm and 𝑑 2 ≈15 mm from the heating wire. The temperature measurements were carried out with two type 𝑁 sheath thermocouples with an outer diameter of 0.5 mm, inserted in grooves with a 0.5 × 0.5mm 2 section machined on the half-length of the sample in order to measure the temperature at the center. The temperature is recorded with 0.1 s time step by a data logger Picolog TC08. The standard deviation of the temperature measurement is 0.01 °C before heating.

The device was placed in an electric fusing kiln (Vecstar furnace) equipped with a temperature controller (Watlow EZ-Zone) ensuring a very good temperature stability inside the kiln (±1 °C).

Two wires were point-welded on the heating wire with a distance 𝑙 w = 235 mm between the two weldings (see Fig. 2). These wires were connected to a voltmeter (Almemo 2890-9) to measure the electrical voltage 𝑈 between the two points. A stabilized power supply (Tektronix PWS2185) produces an electrical current with a constant intensity 𝐼 passing through the wire. These measurements allow the calculation of the electrical resistance per meter of the heating wire with the following formula:

𝑅 el 𝑙 w = 𝑈 𝐼𝑙 w (12)

Estimation method

In previous works [START_REF] Jannot | An improved model for the parallel hot wire: Application to thermal conductivity measurement of low density insulating materials at high temperature[END_REF][START_REF]Refractory Materials-Determination of Thermal Conductivity-Part 2: Hot-Wire Method (Parallel), Standard[END_REF], the Levenberg-Marquardt algorithm [START_REF] Marquardt | An algorithm for least-squares estimation of nonlinear parameters[END_REF] was used to find the values of the unknown parameters 𝜆, 𝜌𝑐 and 𝑅 that minimize the sum of the quadratic errors:

𝑆 = 𝑁 ∑ 𝑛=1 [ 𝑇 exp (𝑑, 𝑡) -𝑇 mod (𝑑, 𝑡) ] 2 (13)
Where 𝑇 mod is calculated by inverse Laplace transform of (2). Using this estimation method, it was found that the estimated values of the thermal contact resistance 𝑅 become negative when the temperature increases which does not make any physical sense. As an example, Fig. 3(a) represents: the experimental curve and the model curve plotted with the triplet (𝜆, 𝜌𝑐, 𝑅) of the estimated values and the residuals (multiplied by 10 for a better reading) for an experiment on the SilPower ® SG at 600 °C for which the estimation over the interval [100 s, 900 s] led to the value 𝑅 = -5.8 K ⋅ W -1 . It can be seen that the estimation residuals are perfectly flat over the estimation interval which means that from 100 s onwards, the Rosseland assumption is verified, hence, in this case, the use of the model is validated. Between 0 s and 100 s (see Fig. 3(b)), the experimental temperature is higher than the model temperature which means that the heat flux transmitted by radiation is higher than that predicted by the model. The negative value of the contact resistance compensates for this model bias and thus has the effect of shifting the temperature so that is obtained 𝑇 exp (100s) = 𝑇 mod (100s). However, it has another undesirable effect illustrated in Fig. 4, which depicts the reduced temperature sensitivity to the contact resistance for different values of the contact resistance. The value 𝑅 = 2 K ⋅ W -1 is an average value consistent with those obtained from COMSOL Multiphysics ® simulations of our device. We can see that the sensitivity to 𝑅 which varies little with the value 𝑅 = 2 K ⋅ W -1 varies significantly for strongly negative values of 𝑅. We will show that this can create an estimation bias on the values of 𝜆 and 𝜌𝑐.

To do so, we performed a simulation based on Eq. ( 2) with the estimated values 𝜆 = 0.22 W m -1 K -1 and 𝜌𝑐 = 6.11 × 10 5 J m -3 K -1 and considering 𝑅 = 2 K W -1 , and then we added to the result a constant value 𝑇 0 . This 𝑇 0 value thus represents the contribution of radiation between 0 and 100 s not taken into account by the Rosseland model.

Two estimation methods were then applied to this curve over the interval [100 s, 800 s]:

• An estimation E1 of (𝜆, 𝜌𝑐, 𝑅) as described above.

• An estimation E2 of (𝜆, 𝜌𝑐, 𝑇 1 ) by calculating the temperature from the (2), considering 𝑅 = 0 K ⋅ W -1 and adding a constant value 𝑇 1 𝑇 0 that is:

𝑇 (𝑑, 𝑡) =  -1 ⎡ ⎢ ⎢ ⎢ ⎣ 𝜑 w 𝑝𝐿 𝐾 0 ( 𝑝 𝑎 𝑑 ) 𝜌 w 𝑐 w 𝜋𝑟 w 2 𝑝𝐾 0 ( 𝑝 𝑎 𝑟 w ) + 2𝜋𝜆 𝑝 𝑎 𝑟 w 𝐾 1 ( 𝑟 𝑝 𝑎 w ) ⎤ ⎥ ⎥ ⎥ ⎦ + 𝑇 1 𝑇 0 ( 14
)
L. Penazzi et al. Where  -1 is the operator of the inverse Laplace transform. 𝑇 0 and 𝑇 1 both stand for the contribution of radiation between 0 and 100s. For the sake of rigor though, these values have been differentiated since 𝑇 0 represents the value added to the preceding equation and 𝑇 1 represents the estimated value of 𝑇 0 with the E2 estimation method. Furthermore, one may state that due to the larger contribution of radiative heat transfer at high temperature, the contact resistance 𝑅 is greatly reduced to approximately zero. Table 1 shows the estimation results for different values of 𝑇 0 . It can be observed that the E2 estimation method leads to a constant error of 0.7 % on 𝜆 and a negligible error on 𝜌𝑐. The E1 method leads to a small estimation error on 𝜆 (less than 1 % if 𝑇 0 < 0.75 °C) but a large error on 𝜌𝑐 that increases with 𝑇 0 . We therefore treated the experimental curves using the E2 estimation method and choosing an estimation interval such that the residuals are flat and centered on zero over this interval.

Validity of the rosseland approximation

To determine the validity limits of Rosseland approximation in the case of the parallel hot-wire method, we performed simulations using the COMSOL Multiphysics ® coupled conduction-radiation transfer solver based on the approximate model P1 [START_REF]Heat Transfer Module User's Guide[END_REF].

Conduction is modeled with the energy balance equation with a radiative source term (Eqs. ( 6) to ( 8)). The radiative source term is modeled with the P1 approximation (Eqs. [START_REF] Zhang | Numerical study of the influence of thermal radiation on measuring semi-transparent thermal insulation material with hot wire method[END_REF] to [START_REF] Marshak | Note on the spherical harmonic method as applied to the milne problem for a sphere[END_REF]). The solution of the transient coupled conduction and radiation equations is solved iteratively by the Finite Element method (FEM). Solid domains formed by the hot wire, the two thermocouples and the sample are meshed with 4974 unstructured tetrahedrals. The maximum size of an element size is 6.36 × 10 -3 m and the minimum size of an element is 3.6 × 10 -5 m.

The mesh is automatically created and adapted for the model's physics settings. The time dependency in the FEM is solved by an implicit BDF (Backward Differentiation Formula) method. This implicit solver uses an backward differentiation formulas with order of accuracy varying from one (also known as the backward Euler method) to two. The time step used for the simulation is 0.1 s.

We first compared the solutions given by COMSOL Multiphysics ® with the exact solution in two asymptotic cases where an analytical solution exists as proposed by Gross and Tran [START_REF] Gross | Radiation effects on transient hot-wire measurements in absorbing and emitting porous media[END_REF]:

(1) The case where the extinction coefficient is very large (𝑘 a = 𝑘 s = 10 8 m -1 ), the temperature evolves as in the case of a purely conductive transfer because 𝜆 r = 0 W ⋅ m -1 ⋅ K -1 .

(2) The case where the material is completely transparent (𝑘 a = 𝑘 s = 0 m -1 ) with perfectly opaque and reflecting surfaces (the heating wire, thermocouple and external border of the sample with 𝜖 = 0) ( non-participant medium).

The geometry considered is a cylinder with an external diameter of 50 mm heated uniformly on its internal diameter of 0.5 mm, both surfaces being opaque and insulated, and the initial uniform temperature being equal to 1000 °C.

As an example, Fig. 5(a) represents the case of the transparent material the curves of the temperature at a distance 𝑑 = 15 mm from the heating wire, obtained with COMSOL Multiphysics ® on the one hand and with the exact analytical model on the other hand for the following values of the parameters: 𝜆 c = 0.5 W ⋅ m -1 ⋅ K -1 , 𝜌𝑐 = 10 6 J ⋅ m -3 ⋅ K -1 . The curves are in very good agreement and estimates made from the curves simulated with COMSOL Multiphysics ® using the ISO 8894-2:2007 method [START_REF] Jannot | Apparent thermal conductivity measurement of anisotropic insulating materials at high temperature by the parallel hot-wire method[END_REF] lead to values of 𝜆 that differ by less than 0.5 % from the nominal value as shown in Fig. 5(b). The ISO 8894-2:2007 method [START_REF] Jannot | Apparent thermal conductivity measurement of anisotropic insulating materials at high temperature by the parallel hot-wire method[END_REF] was detailed in a previous paper [START_REF] Jannot | An improved model for the parallel hot wire: Application to thermal conductivity measurement of low density insulating materials at high temperature[END_REF], it recommends the estimation for values of the ratio 𝑇 (2𝑡)∕𝑇 (𝑡) between 1.5 and 2.4, in the presented case the validity limit of the semi-infinite medium assumption reduces this interval to [1.7, 2.4]. The results obtained in the case where the extinction coefficient is very large are similar.

After this validation phase, we performed simulations with COMSOL Multiphysics ® of the experimental setup shown in Fig. 6. For each simulation, the temperatures of the heating wire and of the two thermocouples placed respectively at 10 mm and 15 mm from the heating wire were recorded.

Preliminary tests were carried out in order to ensure the simulations results independence to temporal and spatial discretization. The simulations were performed between 0 s to 1000 s:

• For two temperatures: 1000 °C and 1500 °C.

• For two typical materials: a low density thermal insulating material (𝜆 c = 0.1 W ⋅ m -1 ⋅ K -1 and 𝜌𝑐 = 2.5 × 10 5 J ⋅ m -3 ⋅ K -1 ) and a high density thermal insulating material (𝜆 c = 0.5 W ⋅ m -1 ⋅ K -1 and 𝜌𝑐 = 10 6 J ⋅ m -3 ⋅ K -1 ).

• For four values of the extinction coefficient 𝛽, ie 500 m -1 , 1000 m -1 , 5000 m -1 and 10 000 m -1 and for each value of 𝛽 the following three values of 𝑘 a were considered: 𝑘 a = 0.0001 × 𝛽; 0.5 × 𝛽; 𝛽 (with

𝑘 s = 𝛽 -𝑘 a )
. The extreme cases correspond respectively to pure scattering and pure absorption.

The absorption coefficient was not chosen equal to zero because it has been shown [START_REF] Jannot | Thermal diffusivity measurement of insulating materials at high temperature with a four-layer (4L) method[END_REF] that in this case the radiative transfer could be represented by a parallel resistance between the two surfaces, which would invalidate the semi-infinite medium hypothesis even at short times. We chose a minimum value equal to 0.0001 × 𝛽.

The refractive index of the medium was chosen equal to unity and the emissivities of the wire and of the thermocouples equal to 0.8. The thermal properties of the heating wire and thermocouple are those of a stainless steel:

𝜆 = 11 W ⋅ m -1 ⋅ K -1 and 𝜌𝑐 = 3.7 × 10 6 J ⋅ m -3 ⋅ K -1 .
The outer surface of the sample was assumed to be opaque, black and insulated.

The temperature curves obtained for the two thermocouples were treated as numerical experiments to estimate the values of the parameters 𝜆 est , 𝜌𝑐 est and 𝑇 1 , using the E2 estimation method described previously. The time interval over which the estimation was performed was adjusted for each curve to obtain flat and zero-centered estimation residuals over this interval. The estimated value 𝜆 est was then compared to the effective conductivity calculated by:

𝜆 = 𝜆 c + 16 3 𝜎𝑇 3 𝛽 ( 15
)
Where 𝜆 c is the phonic conductivity used in the COMSOL Multiphysics ® simulation 𝛽 is the extinction coefficient, and 𝑇 is the average temperature of the measurement point as recommended by Gross and Tran [START_REF] Gross | Radiation effects on transient hot-wire measurements in absorbing and emitting porous media[END_REF]. The estimated value 𝜌𝑐 est was also compared to the nominal value 𝜌𝑐. The analysis of all the results showed that there is no significant difference between the values estimated with 𝑑 1 = 10 mm and those estimated with 𝑑 2 = 15 mm. Fig. 7 shows the relative deviations between the values of effective conductivity 𝜆 est estimated from the temperature at 15 mm and the nominal values calculated by ( 15) as a function of the Stark number defined by Viskanta [START_REF] Viskanta | Heat transfer by conduction and radiation in absorbing and scattering materials[END_REF] and Cintolesi et al. [START_REF] Cintolesi | Numerical simulation of conjugate heat transfer and surface radiative heat transfer using the P 1 thermal radiation model: Parametric study in benchmark cases[END_REF] and presented in ( 16):

𝑁 = 𝛽𝜆 c 4𝜎𝑇 3 (16)
It is proportional to the ratio of the thermal conductivity 𝜆 to the radiative conductivity 𝜆 r . Fig. 7(a) shows that in the case of the PHW method for the estimation error to be less than 2 % it is necessary that:

• 𝑁 > 0.1 in the case of an absorbing and scattering material, see Fig. 7(a).

• 𝑁 > 2 in the case of a purely scattering material, this limit value can be lower for materials of low thermal conductivity (𝑁 > 0.8 for 𝜆 = 0.1 W ⋅ m -1 ⋅ K -1 for example), see Fig. 7(b).

Fig. 8 shows the relative deviations between the values of the volumetric heat capacity estimated from the 15 mm measurement and the nominal values as a function of the Stark number. It shows that for the estimation error to be less than 2 % it is necessary that 𝑁 > 0.5 in the case of an absorbing and scattering material, see Fig. 8(a). In the case of a purely scattering material, the Stark number must be greater than 10 for the accuracy to be better than 5 %.

It should also be noted that the estimation time interval decreases when the extinction coefficient decreases and this is all the more so as the share of scattering is large. On the one hand the time required for radiation to be well taken into account by a conductive model with effective thermal conductivity increases, and on the other hand, the increase of the effective thermal conductivity when the extinction coefficient 𝛽 decreases makes the assumption of the semi-infinite medium less valid. As an example, the estimation interval for the high density thermal insulating material varies from 0 s to 1000 s for 𝛽 = 10000 at 𝑇 = 1000 °C to 80 s to 250 s for 𝛽 = 500 at 𝑇 = 1500 °C.

Results and discussion

Distance calibration

The main source of uncertainty in the estimation of the volumetric heat capacity is the uncertainty in the distance 𝑑 between the heating wire and the thermocouple. The uncertainty on this distance is of the order of 0.5 mm given the brittleness of the materials in which the grooves are made. For a distance 𝑑 = 10 mm, an uncertainty of 0.5 mm leads to uncertainties of respectively 5 % on the distance 𝑑 and 10 % on 𝜌𝑐 [START_REF]Refractory Materials-Determination of Thermal Conductivity-Part 2: Hot-Wire Method (Parallel), Standard[END_REF]. To reduce this uncertainty, we performed three parallel hot-wire experiments at 20 °C for each of the 3 SilPower ® samples which specific heat at 20 °C was measured by differential scanning calorimetry using the Setaram 𝜇dSc 3 apparatus and we obtained: 𝑐 = 760 J -1 kg -1 K -1 . The accuracy of the measurement is estimated at 2 %. These experimental curves were used to estimate the unknown parameters 𝜆, 𝑑 and 𝑅 considering the known value of 𝜌𝑐.

The results are presented in Table 2. A very good reproducibility of the measurements is observed. It is also observed that the deviation between the estimated values and the nominal values (𝑑 nom,1 =10 mm and 𝑑 nom,2 =15 mm) can reach 5 %. Since the uncertainty on 𝜌𝑐 is mainly due to the error on 𝑑, the estimation error on 𝜌𝑐 is reduced to the uncertainty on the reference value used to estimate the distances, i.e. about 2 %. The accuracy of the estimate is thus significantly improved.

Measurements on silpower ®

The estimated values of 𝑑 1 and 𝑑 2 were then used to estimate 𝜆, 𝜌𝑐 and 𝑅 for the experiments performed between 200 °C to 1200 °C. Consistent with the results of the theoretical study, the values of 𝜆 and 𝜌𝑐 estimated with the measured temperatures at distances 𝑑 1 and 𝑑 2 from the heating wire were almost identical. The average deviation observed for 𝜆 is 1.6 % and that for 𝜌𝑐 is 2.1 %. The results obtained for 𝜆 and 𝜌𝑐 are shown in Figs. 9(a) and 9(b) respectively. The points correspond to the average of the 6 measurements performed (3 for each distance) and the bars represent the standard deviations of these groups of 6 points. The 21 standard deviations on 𝜌𝑐 calculated each on 6 measurements were all less than 5 % with an average of 2.6 % which demonstrates the effectiveness of the distance calibration procedure.

The values of thermal conductivity increase as expected with temperature and with density. The values of the specific heat 𝑐 are close for the three materials, which is logical given that they have the same composition. These values are also very close to those measured by Brückner for silica up to 600 °C [START_REF] Brueckner | Properties and structure of vitreous silica. I[END_REF] also shown in Fig. 9. For temperatures above 800 °C, the difference between our values and those of Brückner increases with temperature, these differences will be explained later after the extinction coefficients of the different grades have been estimated. Fig. 10 shows the experimental and model curves from the estimated values as well as the estimation residuals ×10 for the SilPower ® SG at 20 °C and at 1200 °C for the measurement at 𝑑 = 10 mm. It can be The negative value at 20 °C is due to the contact resistance in the absence of significant radiation at this temperature while the positive value at 1200 °C is mainly due to the fact that the Rosseland model underestimates the radiation at short times and that the radiation, negligible at 20 °C, becomes significant at 1200 °C.

N • exp 𝜆 d 1 𝑑 1 -𝑑 nom,1 𝑑 nom,1 𝜆 d 2 𝑑 2 -𝑑 nom,2 𝑑 nom,2 W m -1 K -1 mm % W m -1 K -1 mm % 𝜌 = 581 kg m -3 1

Estimation of the extinction coefficient

The following assumptions are made:

• The thermal conductivities of the solid and air in the porous medium evolve linearly with temperature.

• The radiation is represented by the Rosseland approximation with a constant extinction coefficient 𝛽.

The effective thermal conductivity can therefore be written: For each material, we can then estimate the parameters 𝑎, 𝑏 and 𝛽 by minimizing the sum:

𝜆 = 𝑎 + 𝑏𝑇 + 16𝑛 2 𝜎𝑇 3 3𝛽 (17)
𝑆 = 𝐽 ∑ 𝑗=1 𝐾 ∑ 𝑘=1 [ 𝜆 exp,𝑘 ( 𝑇 𝑗 ) -𝜆 mod ( 𝑇 𝑗 )] 2 (18)
Where 𝜆 mod corresponds to Eq. ( 17), 𝐽 is the number of temperatures for which the measurement of 𝜆 was performed and 𝐾, the number of thermal conductivity estimations carried out at each temperature 𝑇 𝑗 . Six estimations were performed at each temperature (𝐾 = 6), respectively three at a distance 𝑑 1 = 10 mm and three at a distance 𝑑 2 = 15 mm (see Table 2).

The standard deviation of each parameter 𝑎, 𝑏 and 𝛽 is then deduced during the standard procedure of fitting and given by the covariance matrix. The results are presented in Table 3 and the values of the different estimated conductivities are presented in Fig. 11 for SilPower ® SG and SilPower ® HG. The following comments can be made:

• The extinction coefficient decreases very logically when the density decreases.

• The values of the effective conductivity calculated with the estimated values of the parameters 𝑎, 𝑏 and 𝛽 coincide perfectly with the experimental values.

• It is often assumed that the thermal conductivity of the porous solid matrix is proportional to that of the solid [START_REF] Rottmann | Thermal conductivity of evacuated expanded perlite measured with guarded-hot-plate and transient-hot-wire method at temperatures between 295 K and 1073 K[END_REF]. Despite the high uncertainties, the ratios of the thermal conductivity of SilPower ® at 1200 °C to its value at 20 °C (obtained for all three densities) are found to be very close (maximum deviation of 20 %) to that obtained for silica, of which the backbone of SilPower ® is more than 99 % (see Table 3).

L. Penazzi et al. • The accuracy of the 𝛽 estimate decreases when the value of 𝛽 increases: this is explained by the fact that when 𝛽 increases, the share of the radiative conductivity 𝜆 r in the effective conductivity 𝜆 decreases so the sensitivity of 𝜆 to 𝛽 decreases.

kg m -3 m -1 % W m -1 K -1 % W m -1 K -1 % Thermal SiO
Estimates of the extinction coefficients 𝛽 and thermal conductivities 𝜆 c of the solid allow the calculation of the Stark number 𝑁 for each grade and for each temperature. This then allows us to plot the relative differences between the estimated values 𝜌𝑐 and those calculated from the Brückner values as a function of the Stark number 𝑁 (see Fig. 12). We plotted on the same graph the deviations between the values estimated from COMSOL Multiphysics ® simulations for a scattering material (𝑘 s = 0.9999𝛽) and those from Brückner [START_REF] Brueckner | Properties and structure of vitreous silica. I[END_REF]. The deviations calculated from experiments and estimates follow the same trend and are close. SilPower ® thus behaves as a highly scattering material towards radiation which is consistent with its properties: it is a highly reflective material (reflectivity greater than 95 % in the infrared [START_REF]SilPower ® standard and high density[END_REF]) and it contains a volume fraction of 60 % to 70 % of air (transparent) depending on the density. Finally, as an example, a COMSOL Multiphysics ® simulation was performed with the estimated parameters for SilPower ® SG at 1200 °C i.e. 𝛽 = 6920 m -1 ; 𝜆 r = 0.140 W m -1 K -1 ; 𝜆 c = 𝜆 -𝜆 r = 0.263 W m -1 K -1 and considering a highly scattering material (𝑘 s = 0.9999𝛽). This simulation was then treated as an experimental curve to estimate 𝜆 and 𝜌𝑐. The value of 𝜌𝑐 was adjusted to 𝜌𝑐 = 7.00 × 10 5 J m -3 K -1 so that the estimations from the experimental and the simulated curves lead to the same value, 𝜌𝑐 = 8.30 × 10 5 J m -3 K -1 .

Fig. 13(a) shows one experimental curve obtained at 1200 °Celsius, the model curve with the estimated parameters and the estimation residuals. Fig. 13(b) shows the curve simulated with COMSOL Multiphysics ® , the model curve plotted with the estimated parameters and the estimation residuals.

We can see that the residuals have globally the same shape and an amplitude of the same order of magnitude as those represented on Nevertheless, the value of the absorption coefficient used to match the COMSOL Multiphysics ® simulation and the experiment is not realistic (𝑘 a = 6.9 m -1 ). This may be due to the approximation of the P1 model or to the hypothesis of a gray medium. It will therefore be necessary to develop a more accurate model if we want to represent the whole experimental curve in order to estimate for example the extinction or diffusion coefficient of the medium.

Conclusion

The influence of the radiation heat transfer on the estimation of the thermal conductivity of a semi-transparent material by the parallel hot-wire method has been first investigated with COMSOL Multiphysics ® using the coupled conduction-radiation model based on the P1 approximation for the radiation calculation. It was shown that the thermal conductivity may be estimated with an accuracy better than 2 % if the Stark number 𝑁 (representing the ratio conduction/radiation) is:

• Greater than 0.1 for a purely absorbing or absorbing and scattering medium.

• Greater than 2 for a purely scattering medium.

It has also been shown that the classical estimation method based on the estimation of the three parameters (𝜆, 𝜌𝑐) leads to negative values of the thermal contact resistance 𝑅 ( which makes no physical sense) when radiation increases and creates an estimation bias. An new estimation method has been proposed to avoid this estimation bias.

A calibration process of the distance between the hot wire and the thermocouple has also been proposed and validated, enabling a more accurate estimation of the volume heat capacity. An experimental study carried out on three insulating materials with densities ranging from 581 kg m -3 to 910 kg m -3 and at temperatures ranging from 20 °C to 1200 °C confirms the results of the theoretical study. A method



  2022.107690 Received 4 January 2022; Received in revised form 2 May 2022; Accepted 6 May 2022 Nomenclature Greek Symbols 𝛽 Extinction coefficient (m -1 ) 𝜆 Thermal conductivity of the considered material (W m -1 K -1 ) 𝛺 Solid angle (sr) 𝜙 Heat flux density (W m -2 ) 𝜌 Density of the material (kg m -3 ) 𝜎 Stefan-Boltzmann constant (W m -2 K -4 ) 𝜏 Optical thickness (-) 𝜃 Laplace transform of 𝑇 (𝑑, Radiant intensity (W m -2 sr -1 ) 𝒏 Outward unit direction on the surface (-) 𝑎 Thermal diffusivity of the isotropic material (m 2 s -1 ) 𝑐 Heat capacity of the material (J kg -1 K -1 ) 𝑑 Distance between the heating wire and the thermocouple (m) 𝐺 Incident radiation (W m -2 ) 𝐼 Electrical intensity passing through the heating wire (A) 𝑘 a Absorption coefficient (m -1 ) 𝑘 s Scattering coefficient (m -1 ) 𝐾 𝑖 Modified Bessel functions of the second kind 𝑖 ∈ {0; 1} 𝑙 w Length of the heating wire (m) 𝑁 Stark number (-) 𝑛 Refractive index (-) 𝑃 Phase function (-) 𝑝 Laplace parameter (s -1 ) 𝑄 r Radiative heat source (W m -3 ) 𝑅 Thermal contact resistance between the heating wire and the sample (K W -1 ) 𝑟 w Radius of the heating wire (m) 𝑅 el Electrical resistance of the heating wire (Ω)

Fig. 1 .

 1 Fig. 1. Schematic view of a parallel hot-wire set-up.

Fig. 2 .

 2 Fig. 2. Schematic view of the measurement set-up.

Fig. 3 .

 3 Fig. 3. Experimental and simulated curves with estimation residuals ×10 for an experiment with SilPower ® SG at 600 °C.

Fig. 4 .

 4 Fig. 4. Reduced sensitivity of the temperature according to thermal contact resistance 𝑅 values (𝜆 = 0.22 W m -1 K -1 ; 𝜌𝑐 = 6.11 × 10 5 J m -3 ⋅ K -1 ).

Table 1 11 𝑇 1 °C

 1111 Estimated values of the parameters by the two methods E1 and E2 for different values of the temperature step due to radiation between 0 s to 100 s with 𝜆 = 0.22 W m -1 K -1 , 𝜌𝑐 = 6.11 × 10 5 J m -3 K -1 . 1 K -1 0.220 0.219 0.219 0.218 0.217 0.217 0.216 𝜌𝑐 (×10 5 ) J m -3 K -1 60 -9.8 -15.4 -20.8 -26.0 -31.3 E2 𝜆 W m -1 K -1 0.219 𝜌𝑐 (×10 5 ) J m -3 K -1 6.

Fig. 5 .

 5 Fig. 5. (a) -Temperature curves simulated with COMSOL Multiphysics ® and with the exact analytical model (𝜆 c = 0.5 W m -1 K -1 ; 𝑘 a = 𝑘 s = 0 m -1 ) ; (b) -Estimated values of 𝜆 by processing the COMSOL Multiphysics ® curve with the standard method ISO 8894-2:2007.

Fig. 6 .

 6 Fig. 6. Diagram of the modeled system with 𝑑 1 = 10 mm and 𝑑 2 = 15 mm.

Fig. 7 .

 7 Fig. 7. Relative deviation ( 𝛥𝜆 𝜆 = 𝜆 est -𝜆 nom 𝜆 nom ) between the effective conductivity values (𝜆 est ) estimated from the 15 mm PHW measurements and the nominal values (𝜆 nom = 𝜆 c + 𝜆 r ) as a function of the Stark number in (a) the absorption-diffusion case and (b) the pure scattering case.

Fig. 8 .

 8 Fig. 8. Relative deviation ( 𝛥𝜌𝑐 𝜌𝑐 = 𝜌𝑐 est -𝜌𝑐 0 𝜌𝑐 0 ) between the volumetric heat capacities values estimated (𝜌𝑐 est ) from the 15 mm PHW measurements and the nominal values (𝜌𝑐 0 ) as a function of the Stark number in (a) the absorption-diffusion case and (b) the pure scattering case.

Fig. 9 .

 9 Fig. 9. (a) -Estimated values of the effective thermal conductivity ; (b) -Estimated values of the specific heat and values given by Brückner [22]. (SG: Standard density, MG: Medium density, HG: High density).

Fig. 10 .

 10 Fig. 10. Experimental curve and simulated curve with the values estimated between 100 s and 800 s both with estimation residuals ×10 for SilPower ® SG (a) at 20 °C and (b) at 1200 °C for 𝑑 = 10 mm.

Fig. 11 .

 11 Fig. 11. Estimated values of the various conductivities for SilPower ® SG (Standard density) and SilPower ® HG (High density).

Fig. 12 .

 12 Fig. 12. Relative deviation between the values of the specific heat estimated from experiments or simulations and those given by Brückner [22].

Fig. 13 .

 13 Fig. 13. (a) -Experimental curve and model curve ; (b) -Simulated curve with COMSOL Multiphysics ® and model curve.

Fig. 9 ( 3 = 15 .

 9315 b), which confirms the very scattering character of SilPower ® . Moreover, the deviation between the nominal and estimated values of 𝜌𝑐 which is 𝛥𝜌𝑐 𝜌𝑐 = 8.3-7.0 8.6% is of the same order of magnitude as the one observed with the experiments.

Table 2

 2 Estimated values of SilPower ® thermal conductivity and heating wire-thermocouple distances at 20 °C.

  which validates the model of an effective thermal conductivity on this time domain. The estimation at 20 °C leads to perfectly flat residuals as early as 40 s while at 1200 °C we observe a divergence of the residuals up to 100 s. The estimated values of the temperature shift are 𝑇 1 = -0.039 °C at 20 °C and 𝑇 1 = 1.02 °C at 1200 °C.

		2 3	0.123 0.124 0.125	9.99 10.00 10.01		0.121 0.124 0.125	14.65 14.67 14.67
		Mean 0.124	10.00 0.0	0.123	14.66 2.27
	𝜌 = 741 kg m -3	1 2 3	0.169 0.169 0.168	9.72 9.72 9.71		0.160 0.160 0.158	14.17 14.18 14.78
		Mean 0.169	9.72	2.80	0.159	14.17 5.53
	𝜌 = 910 kg m -3	1 2 3	0.213 0.217 0.214	9.38 9.43 9.46		0.211 0.217 0.214	14.78 14.91 14.88
		Mean 0.215	9.42	5.80	0.214	14.85 0.96
	seen that the residuals are perfectly flat and centered on the estimation interval [100 s; 800 s]

Table 3

 3 Estimated values of parameters 𝑎, 𝑏 and 𝛽 and their uncertainties.

	𝜌	𝛽	𝜎 𝛽	𝑎	𝜎 𝑎	𝑏	𝜎 𝑏	𝜆 𝑐 1200 °C 𝜆 𝑐 20 °C

Data availability

Data will be made available on request.

L. Penazzi et al. enabling the estimation of the extinction coefficient from thermal conductivity measurements at various temperatures has been presented and successfully applied to the three tested materials.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.


Negin Alisoltani
email: negin.alisoltani@geotwin.io

Mostafa Ameli
email: mostafa.ameli@univ-eiffel.fr

Mahdi Zargayouna
email: mahdi.zargayouna@univ-eiffel.fr

Ludovic Leclercq
email: ludovic.leclercq@univ-eiffel.fr

A SHAREABILITY CLUSTERING METHOD TO SOLVE THE DYNAMIC RIDE-SHARING PROBLEM CONSIDERING NETWORK CONGESTION

Keywords: Dynamic ride-sharing, Clustering, Shareability function, Optimal fleet management

One of the most challenging facets of dynamic ride-sharing problems is the automated matching process, especially in large-scale networks where the number of received requests is huge. Matching users to trips is very challenging in real-time since it must happen very quickly while it must ensure high quality solutions. In this paper, we rely on an exact solving method to ensure the quality of the solution, while using AI-based techniques to limit the number of requests that we feed to the solver. More precisely, we propose a clustering method based on a new space-time shareability function to put the most shareable trips inside separate clusters. Each cluster is then managed with a proposed heuristic framework in order to solve the matching problem inside each cluster. To validate our proposal, we employ the proposed method on the network of Lyon city in France, with half-million requests in the morning peak from 6 to 10 AM. The results demonstrate that the algorithm can provide high-quality solutions in a short time for large-scale problems.

INTRODUCTION

Ridesharing originated as a general concept where individual travelers share a vehicle for a trip and split travel costs with others that have similar itineraries and time schedules (Furuhata et al.,[START_REF] Furuhata | Ridesharing: The state-of-the-art and future directions[END_REF]. The spread of mobile devices and the development of Global Positioning System (GPS) make it possible for all the transport operators to adapt in real-time the transportation supply to travel demand (Zargayouna et al.,[START_REF] Zargayouna | Multiagent simulation of realtime passenger information on transit networks[END_REF]. These options make it possible to have access to the vehicles' position at any time and perform the matching process of ride-sharing in real-time. These possibilities have led to the development and progress of a new type of ride-sharing which is called real-time ride-sharing, also known as dynamic ride-sharing or ad hoc ride-sharing. Dynamic ride-sharing refers to a system that supports an automatic ride-matching process between participants on very short notice or even en-route (Agatz et al.,[START_REF] Agatz | Dynamic ride-sharing: A simulation study in metro Atlanta[END_REF].

The first issue in this application is matching the passengers and the fleet vehicles by solving a fleet management problem. Effective and efficient optimization technology that matches drivers and riders in real-time is one of the necessary components for a successful dynamic rideshare system (Agatz et al.,[START_REF] Agatz | Optimization for dynamic ride-sharing: A review[END_REF]. Matching users to trips is very challenging in real-time since it must happen very quickly. Clustering methods enhance the scalability of matching methods by reducing the search space, making it possible to parallelize and balance the computing workload, and as a consequence speeding up the matching computations.

In this paper, we propose a clustering method based on a "Shareability Function" (SF) that considers all the trips' possible matching situations. Two trips can be shared either in parallel (passengers are in the same vehicle at some point) or in sequence (Alisoltani et al.,[START_REF] Alisoltani | A Sequential Clustering Method for the Taxi-Dispatching Problem Considering Traffic Dynamics[END_REF]. The shareability function computes the extra travel time that the vehicle has to spend to service each matching situation, compared to the situation where each trip is serviced independently without sharing. Then the shareability function for each two trips is the minimum computed value among the three different situations (parallel, sequence and independence). So, when the shareability function is low for two specific requests, it means that they have high potential to be shared efficiently. We propose a clustering method based on this new function to put the most shareable trips in separate clusters. We propose a heuristic method to solve the matching problem inside each cluster. The final algorithm can provide high quality solutions in a short time for large-scale problems.

It is crucial to consider the interactions between the system and the network traffic in dynamic ride-sharing system evaluation. Also, it is important to consider other vehicles like personal cars besides the ride-sharing service vehicles in the network. In this paper, we define two different models to simulate the functioning of the proposed dynamic ride-sharing system: The "plant model", based on Macroscopic Fundamental Diagram (MFD), is used to simulate the real traffic conditions and considers both service vehicles and personal vehicles in the network; The "prediction model", based on the current mean speed, is used to calculate the travel times during the assignment process (Alisoltani et al.,Alisoltani et al.,[START_REF] Alisoltani | A Sequential Clustering Method for the Taxi-Dispatching Problem Considering Traffic Dynamics[END_REF][START_REF] Alisoltani | Optimal fleet management for real-time ride-sharing service considering network congestion[END_REF].

To validate our approach, we employ the proposed method on the real data of Lyon city in France with half-million requests in the morning peak from 6 to 10 AM. To demonstrate the clustering method's performance, we have compared the proposed method with two other clusteringbased methods from the literature (spatial and temporal clustering). The results show that the method proposed in this paper can make a significant improvement in the quality of the solution and the computation time.

Since ride-shares are established on-demand, a ride-sharing system is similar to other ondemand forms of passenger transit such as taxis and dial-a-ride services like airport shuttles (Fu-TRB 2022 Annual Meeting Original paper ruhata et al., 1). The clustering method proposed in this study can be used for other similar systems.

The remainder of this paper is structured as follows. Next section is the literature review on ride-sharing and clustering methods. In the third section, we explain our clustering method for dynamic ride-sharing. Then we define our travelers matching algorithm and describe the two simulation models used as plant and prediction models respectively. I section "Experiments" we present the case study we used and provide the results of our experiments. The last section concludes this paper.

LITERATURE REVIEW

In this section, we review the latest studies on different solution methods and clustering approaches for real-time and dynamic ride-sharing.

Dynamic ride-sharing problem addresses short-term matching or even en-route matching (Agatz et al.,[START_REF] Agatz | Dynamic ride-sharing: A simulation study in metro Atlanta[END_REF]. This fact makes the assignment problem more complex. In some studies, researchers try to narrow the feasible solution space to make the computations faster and be able to assign the vehicles to the requests that are coming at each time to the system. For example, Liu et al. [START_REF] Liu | A branch-and-cut algorithm for a realistic dial-a-ride problem[END_REF] present a method to tighten travelers time windows and eliminate unnecessary variables and constraints to narrow the solution space.

Recently Mourad et al. [START_REF] Mourad | A survey of models and algorithms for optimizing shared mobility[END_REF] has presented a survey of models and algorithms for optimizing shared mobility, and they have shown that one of the most critical problems in the solution for these systems is computation time and the solutions quality.

Compared to these approaches, this paper proposes a method that combines the benefits of exact solving in terms of solution quality and the benefits of heuristics and meta-heuristics in terms of computation time and scalability. On the one side, we present an algorithm based on branch-and-bound that can provide the optimal solution for small instances of the problem. On the other side, we introduce a method based on clustering that groups the requests and put the most shareable trips in the same clusters. The exact algorithm is executed with the requests inside each cluster. A rolling horizon approach is introduced to handle the trip requests in real-time. We show that the clustering method provides high-quality solutions while reducing the computation time significantly.

The number of request for the mobility service at every time is huge in the large-scale problems and the complexity of the solution method for the matching problem increases exponentially by small increase in the service demand. [START_REF] Qiang | Clustering algorithm for urban taxi carpooling vehicle based on data field energy[END_REF] propose an algorithm that uses the dataset of taxi get-off points and performs a clustering of taxis on urban roads. They compare their method with classical clustering methods.

However, the taxi clustering data in their study are conducted in a static environment. In Chen et al. [START_REF] Chen | Solving the first-mile ridesharing problem using autonomous vehicles[END_REF], all pickup points are partitioned into several clusters, the vehicles dispatching and the ride-TRB 2022 Annual Meeting Original paper sharing problem are solved in each cluster. Bard and Jarrah (16) show that an appropriate solution for large-scale problems is clustering the demand nodes and downsizing the network. To speed up the computation, Some researches try to limit the feasible region with clustering methods. In these methods, the demand nodes in the network are usually divided into geographically dense clusters (Özdamar and Demir,Sáez et al.,[START_REF] Özdamar | A hierarchical clustering and routing procedure for large scale disaster relief logistics planning[END_REF][START_REF] Sáez | Hybrid adaptive predictive control for the multivehicle dynamic pick-up and delivery problem based on genetic algorithms and fuzzy clustering[END_REF]. Li and Chung [START_REF] Li | Ride-sharing under travel time uncertainty: Robust optimization and clustering approaches[END_REF] propose an extended insertion algorithm in conjunction with a Tabu search method, and a cluster-first-route-second approach to find heuristic solutions. Under the assumption of one person per ride, Alonso-Mora et al. [START_REF] Alonso-Mora | On-demand highcapacity ride-sharing via dynamic trip-vehicle assignment[END_REF] show that 98% of the taxi rides currently served by over 13,000 taxis could be served with just 3,000 taxis of capacity four in New York City. The proposed clustering method in this paper, seeks the potential trips to be shared in the same vehicle and makes sure that the matching algorithm can use the maximum capacity of a vehicle to serve the requests for a specific route.

In a recent study, Yu et al. [START_REF] Yu | Optimal routing of multimodal mobility systems with ride-sharing[END_REF] consider the combination of the ride-sharing and public transportation services and formulate a mixed integer programming model for the multimodal transportation planning problem. They propose a heuristic approach (i.e., angle-based clustering algorithm) and compare its efficiency with the exact solution for different settings. They show that the clustering algorithm works well both in small and large settings. However, it is crucial to consider the interactions between the system and the network traffic in the ride-sharing problem.

In this paper, we define two distinct models for the dynamic ride-sharing to calculate the travel times during the assignment and simulate real traffic conditions.

CLUSTERING METHOD FOR DYNAMIC RIDE-SHARING Shareability function

To perform the clustering on the requests (N t ) received by the system at time t, we define the "Shareability Function" (SF i, j ) between request i and request j (∀i, j ∈ N t ). This function computes the difference between the travel time when the two trips are shared and the travel time to serve each trip individually, for each pair of requests. In case of sharing, we consider three different situations for each pair of trips (Fig. 1). In the first situation (1), two trips can be shared, and the vehicle first drops off the first passenger, and then it goes to the second passenger drop off (destination) point. In this situation, the travel time for the first passenger (T p i ) would be the summation of her/his waiting time (W s i ), the travel time between her/his origin and the next passenger origin (T OO i, j ) and the travel time between the next passenger origin and her/his destination (T OD j,i ). Similarly, the travel time for the second passenger (T p j ) would be the summation of her/his waiting time (W s j ), the travel time between her/his origin to the first passenger's destination (T OD j,i ) and the travel time from the first passenger's destination and her/his destination (T DD i, j ). So, SF i, j (shareability function for trips i and j) when two trips have the first situation can be computed as in equation 1.

T p i + T p j = W s i + T OO i, j +W s j + T OD j,i + T DD i, j ∀i, j ∈ N SF 1 i, j = T p i + T p j -(T OD i,i + T OD j, j +W i +W j ) (1)
In situation (2), the second passenger is served by the service vehicle while the first passenger is on board. Thus, the travel time for the second passenger (T p j ) is the same as when served individually (summation of the waiting time and the time from the origin to the destination) and the travel time for the first passenger is the travel time of all the links from the first stop point (O i )

to the last one (D j ):

T p i + T p j = W s i + T OO i, j +W s j + T OD j, j + T DD j,i ∀i, j ∈ N SF 2 i, j = T p i + T p j -(T OD i,i + T OD j, j +W i +W j ) (2)
In the third situation (3) that we consider for two trips, the trips are not shared, but the vehicle can serve two passengers sequentially (Alisoltani et al.,[START_REF] Alisoltani | A Sequential Clustering Method for the Taxi-Dispatching Problem Considering Traffic Dynamics[END_REF]. This situation must be considered in the shareability index in order to encourage putting these trips in the same group while solving the optimization problem. The travel time for both passengers in this situation is the same as when they are served individually. But the vehicle travel time can decrease if the travel time between the first destination and the second origin is less than the summation of the travel time between the first origin and the closest depot and the travel time between the start depot and the second origin.

T p i + T p j = W s i + T OD i,i +W s j + T OD j, j ∀i, j ∈ N SF 3 i, j = T p i + T p j -(T OD i,i + T OD j, j +W i +W j ) SF 3 i, j = W s i +W s j -(W i +W j ) (3)
The best situation for each two trips is the situation with minimum SF. So the algorithm chooses the condition that the additional travel time is minimum for sharing each pair of trips:

SF i, j = min{SF 1 i, j , SF 2 i, j , SF 3 i, j } (4)
TRB 2022 Annual Meeting Original paper

Clustering based on a dissimilarity function

After computing the shareability function, we have the function value for each pair of requests that creates a "shareability matrix". The shareability matrix is a dissimilarity matrix for the received requests: the higher the value, the least likely the two requests will be served together. In our approach, this matrix is used in the clustering process. We perform the clustering using the computed dissimilarity matrix. When we make clusters based on the SF, we put in the same cluster the trip requests that have more potential to be shared (the trips that have lower SF). There are two main categories of algorithms for clustering based on the shareability matrix, both are potentially relevant for our large-scale ride-sharing application: In Qiang and Shuang-Shuang ( 14), the authors have compared the hierarchical clustering and kmeans clustering for urban taxi carpooling in a static environment. They show that compactness and separation are almost the same for hierarchical and k-means clustering for large cluster sizes.

However, in dynamic large-scale problems, the results might be different. Besides, computation time becomes the critical criterion in this context. Thus, in our experiments, we implement both clustering methods and we assess their performance, considering both solutions quality and computation time.

Multidimentional scaling and k-means clustering

K-Means method is a partitional clustering approach for decomposing the problem into independent subsets. It defines clusters of data based on their similarity, minimizing within-cluster variances. In the clustering procedure, the preferred number of clusters (k) should be specified in the algorithm before execution. The most common algorithm uses an iterative refinement technique.

A common initialization step associates each observation (data point) with a cluster, and computes the set of cluster centroids as the mean of the observations of the cluster. Afterward, iterations will serve to optimize the clusters.

k-means clustering takes place based on the distance between points. Based on the study in Paea and Baird [START_REF] Paea | Information Architecture (IA): Using Multidimensional Scaling (MDS) and K-Means Clustering Algorithm for Analysis of Card Sorting Data[END_REF], to be able to apply the appropriate k-mean clustering method, we need to convert the dissimilarity matrix into a distance matrix. Therefore, we use the multidimensional scaling method and change the shareability matrix into a distance matrix (Wang and Boyer,[START_REF] Wang | Feature learning by multidimensional scaling and its applications in object recognition[END_REF].

Note that we have implemented the method using the mathtoolbox in C++).

After extracting the distance matrix, we are able to create the same size clusters for the data received at every assignment time step and put the most shareable trips in separate clusters using the modified k-means clustering method. This k-mean clustering is used to obtain clusters in preferred sizes (Ganganath et al.,[START_REF] Ganganath | Data clustering with cluster size constraints using a modified k-means algorithm[END_REF]. Accordingly, considering the objective function value, we can find the best trade-off between cluster size and computation time. Putting the requests in clusters can also provide the opportunity for parallel computations. So the problem is divided into TRB 2022 Annual Meeting Original paper multiple sub-problems, and we favor a uniform distribution of requests among clusters to decrease computation times and facilitate parallel computations of each sub-problem.

Hierarchical clustering

Hierarchical clustering offers a flexible and no-parametric approach and is an algorithm that builds hierarchy of clusters (Murtagh,[START_REF] Murtagh | A survey of recent advances in hierarchical clustering algorithms[END_REF].

We use the agglomerative hierarchical method, which starts with taking singleton clusters (that contain only one request per cluster) at the bottom level and continue merging two clusters at a time to build a bottom-up hierarchy of the clusters (Reddy and Vinzamuri, 30).

MATCHING ALGORITHM FOR DYNAMIC RIDE-SHARING

The matching algorithm aims to identify the travelers who can share their trips and assign them to a vehicle. In the assignment, it is important to consider both passengers' and service provider's objectives. According to the state-of-the-art, the most important operation objective for the service provider is to minimize the total travel time, and the total travel distance of vehicles (Ota et al., Qian et al.,[START_REF] Ota | STaRS: Simulating taxi ride sharing at scale[END_REF][START_REF] Qian | Optimal assignment and incentive design in the taxi group ride problem[END_REF]. The passengers also need to get to the destination on time and have the minimum waiting time (Hyland and Mahmassani,[START_REF] Hyland | Sharing Is Caring: Dynamic Autonomous Vehicle Fleet Operations Under Demand Surges[END_REF]. So we define the objective function for the matching algorithm as below where i is the index of passengers, and m is the index of vehicles, W i is the waiting time for passenger i, T p i is the travel time for passenger i, T v m is the travel time for vehicle m, D m is the travel distance for vehicle m and α, β , γ, and δ are the weights of each objective after normalization.

min ∑ i∈P (α.W s i + β .T p i ) + ∑ m∈M (γ.T v m + δ .D m ) (5)
The main constraints for the matching process are capacity constraints, time constraints, and assignment constraints. In addition, we have considered limitations on the number of sharing.

The number of sharing for each passenger defines the allowed number of other passengers that can share their trip with this current passenger. So when the number of sharing is zero, it means that the service is not able to share and we just consider the in-sequence trips in the clustering function. But when we increase the number of sharing to one, it means that the service can serve 2 passengers at the same time. Also, we made sure that there is a sufficient number of vehicles in the fleet to serve all the requests.

The vehicle can be in two different situations. It is either circulating in the network to serve the on-board passengers (en-route vehicle) or waiting at a stop location for the new passengers (idle vehicle). Our proposed algorithm considers both situations. When the system receives a new request, it first checks the en-route vehicles that have not yet reached their maximal occupancy.

So, the first part of the algorithm assigns the new requests in priority to en-route vehicles. The algorithm works to minimize the total travel time for both vehicles and passengers while respecting the capacity of the vehicle's constraint and passengers' time window constraint. For the remaining requests, the second part of the algorithm checks the idle vehicles waiting in depots to assign the passengers. For this step, we introduce an exact method based on the branch-and-cut concept.

The method creates branches of routes and then chooses the optimal solution among the feasible routes. Fig. 2 shows this part of the algorithm. Then we introduce a rolling horizon method to solve the problem dynamically. The requests that are assumed to be known are the only ones over the next rolling horizon (20 minutes). We also use an "assignment time horizon" of 10 minutes.

TRB 2022 Annual Meeting Original paper

That means that every 10 minutes, we execute our optimization algorithm while considering the requests of the next 20 minutes. We present a small example with four requests to show the execution of the algorithm. Table 1 presents the requests. Each request has an associated number of requested seats, a maximum number of accepted sharing, an earliest pick-up time and a latest arrival time. The algorithm starts by creating branches of routes to serve the requests. It sends a car from the closest stop location to pick the passengers up at the origin point, and then it continues by adding the feasible points to the branches. First, the algorithm finds feasible branches. Figure 3 shows the final feasible solutions for the problem. The algorithm can find four feasible solutions for the problem. Table 2 shows the total passengers' waiting time, total vehicles' travel distance, and the number of vehicles for these solutions (note that the waiting time for passengers is the difference between the passengers' pick up time and the passengers' desired departure time). Solution A is when a car serves each passenger separately without sharing. In this solution, the waiting time is minimum, and the passengers wait for 28 seconds to be picked up. Solution B serves passengers The ride-sharing service's optimization system uses estimates for the predicted travel time obtained from a so-called "prediction model". When the fleet management plan is executed, a gap usually exists between the estimation and the real traffic condition. The so-called "plant model" requires dynamic adjustment of the initial assignment to fit with the conditions observed, and we use it to represent the real traffic condition.

To provide a realistic service, in the proposed simulation component of the dynamic ridesharing system, we accurately distinguish the prediction and the plant models.

The trip-based MFD is used as the plant model (Ameli et al.,[START_REF] Ameli | Departure Time Choice Models in Congested Transportation Systems Based on Mean Field Games[END_REF] to consider individual trips while keeping a very simple description of traffic dynamics (Lamotte and Geroliminis, Ameli,[START_REF] Lamotte | The morning commute in urban areas: Insights from theory and simulation[END_REF][START_REF] Ameli | Heuristic Methods for Calculating Dynamic Traffic Assignment[END_REF]. The general principle is to derive the inflow and outflow curves. When n(t) is the number of en-route vehicles at time t and the mean speed of travelers is V (n(t)) at every time t, the travel distance L i by a car i entering at time t -T (t) must satisfy the following equation:

L i = t t-T (t) V (n(s))ds (6)
The function V (n(t)) is the speed macroscopic fundamental diagram and can be derived from common observations for a transportation network (Leclercq et al.,[START_REF] Leclercq | Macroscopic fundamental diagrams: A cross-comparison of estimation methods[END_REF]. For more details on the functioning of trip-based MFD, readers can refer to (Mariotte and Leclercq,[START_REF] Mariotte | Flow exchanges in multi-reservoir systems with spillbacks[END_REF].

The prediction model estimates the traffic situation for the next assignment time horizon (every 10 minutes) to carry out travel time prediction during optimization. In our method, we assign the passengers to the cars based on this prediction. The prediction model is based on the direct travel time from each point i to j based on the current mean speed and the associated shortest path between the two points for the next 10 minutes. Then the optimization algorithm assigns all the requests for the next 10 minutes to the en-route cars or empty waiting cars.

In the rolling horizon method, the assignment procedure rolls over a specific horizon for the requests announced of a particular optimization step. As stated earlier, the rolling horizon in TRB 2022 Annual Meeting Original paper this paper is 20 minutes, and the optimization time step is 10 minutes. So, the requests of the next 20 minutes that have not yet been assigned are optimized every 10 minutes. Some requests are re-optimized every 10 minutes. If a trip has been assigned to a vehicle which has left the depot, the algorithm does not assign it again, but if a trip is in the schedule of a waiting vehicle, the algorithm places it in the set of optimized trips in the particular horizon and re-optimizes it.

In this method, every T H 2 time step, we stop the simulations and solve a new assignment problem for the requests over a new full rolling horizon (T H). Some requests may arrive just after the end of a simulation period, and this method prevents the system from being myopic to the new demand.

EXPERIMENTS Case study

In this paper, we use data from the network of Lyon city (the second-largest urban area of France)

to evaluate the impact of proposed clustering method. The network area is more than 80 km 2 and the origins/destinations (ODs) set contains 11,314 points (Ameli et al., 39). Fig. 4 shows the network of the city.

FIGURE 4 : Network of Lyon

The network is loaded with travelers of all ODs with a given departure time to represent the morning peak hour from 6 AM to 10 AM. The number of trips during this period is 484,690. We have 279,382 personal trips in the network and 205,308 demand for the service cars in the system.

The service provider has a fleet of vehicles in a ride-sharing system to serve the service requests. Participating service vehicles start up from a number of known locations or depots and after serving the assigned requests, they stop at this allowed locations to wait for the next passengers.

TRB 2022 Annual Meeting Original paper

Determining the proper clustering method

As we said in section 4.2, we can use either the k-means clustering method or the hierarchical method to cluster the requests based on the presented shareability matrix. In our method, both the quality of the clustering method and the computation time are very important. The time complexity of k-means is linear, while that of hierarchical clustering is quadratic. Besides, k-means clustering requires prior knowledge of number k of clusters and also needs to use mutidimensional scaling to convert the similarity matrix into the distance matrix. On the other side, we can stop at whatever number of clusters we find appropriate in hierarchical clustering by interpreting the dendrogram. As we use the agglomerative hierarchical method, we can have larger clusters faster with the hierarchical method. To choose the best clustering method, we have compared both methods considering the quality of objective function and the computation time for different sizes of problems.

Table 3 shows the objective function and computation time with k-means clustering and hierarchical clustering method for different sizes of problems.

To have a baseline for the comparisons, we have computed the optimal objective function (solution method without clustering) for 112 and

The size of clusters

We try to have the same size clusters (to avoid too big or too small clusters) to keep the computation time low and have the opportunity for sharing in all the clusters.

There are different methods in the literature to choose the optimal size of clusters. Here, the quality of the clusters (how similar are points within a cluster) is very important. Furthermore, we have to be sure that the clusters are separated from each others, and the possibility of sharing two trips from two different clusters is minimum. Thus the best way to find the optimal size of clusters is to use the Sum of Squares method (SS) (Krzanowski and Lai,[START_REF] Krzanowski | A criterion for determining the number of groups in a data set using sum-of-squares clustering[END_REF]. It is a clustering validation method that chooses the optimal size of clusters by minimizing the Within-cluster Sum of Squares (WSS) (a measure of how tight each cluster is) and maximizing the Between-cluster Sum of Squares (BSS) (a measure of how separated each cluster is from the others). We compute the WSS and BSS for all the clusters in different periods to evaluate the optimal size of clusters in different demand situations.

K-means clustering

We compute the sum of squares for k-means clustering with cluster sizes from 10 to 50. Fig. 5 shows that when the size of clusters is 30, we can find the best trade-off between WSS and BSS.

So, we choose 30 as the size of clusters when we want to cluster the service requests with k-means clustering.

Hierarchical clustering

In the hierarchical method, we can keep the first cluster of the desired size at the bottom of the dendrogram to have the same size clusters (Hippocamplus, 41). The crucial point here is to find an approximation for the size of clusters in the hierarchical method. As we have explained earlier, we use the hierarchical clustering when the number of demand is huge. So, we analyze the size of clusters from 75 to 300. Also, in large-scale networks with high levels of traffic congestion, the demand density is different in different times of a day (during the on-set and off-set of congestion).

TRB 2022 Annual Meeting Original paper Fig. 6 shows the SS method at different times of the simulations. Increasing the size of clusters decreases the BSS. It means that more number of clusters can ensure that the clusters are separate from each others. We have determined the cluster sizes that minimize the WSS and maximize the BSS. At 6 and 7 AM, the cluster sizes 100 and 125 can make this trade-off between WSS and BSS. At 8 and 9 AM, the best cluster sizes are 125 and 150. Therefore, the best cluster size to do the simulations for this scale is 125.

Comparing the shareability clustering method with spatial and temporal methods Some researches on ride-sharing use clustering methods to scale-down the matching problem and compute the assignment for the service vehicles faster. They usually divide the space geographically and use a spatial clustering to downsize the problem (Chen et al.,[START_REF] Chen | Solving the first-mile ridesharing problem using autonomous vehicles[END_REF]. The important factor in spatial clustering is the distance between the trips' origins. So two corresponding trips can be in the same cluster if the distance between their origins is small. Another approach is to cluster the trips based on the time in a temporal clustering method. In the temporal clustering, we put two TRB 2022 Annual Meeting Original paper

FIGURE 1 :

 1 FIGURE 1 : Trip situations

1 .

 1 Partitional clustering algorithms (Celebi, 24) cluster the data into k clusters. One of the usual algorithms for partitioned clustering is k-means clustering. K-means clustering is simple, fast, and flexible. 2. Hierarchical clustering methods in which the clusters are arranged in a tree-like structure. Hierarchical clustering can be divided into Agglomerative hierarchical clustering (AHC) and divisive clustering (Vora et al., 25).

1 and 4 FIGURE 2 :

 42 FIGURE 2 : Final assignment algorithm

FIGURE 3 :

 3 FIGURE 3 : Feasible routes for the example

FIGURE 5 :

 5 FIGURE 5 : Sum of squares method for finding the optimal size of clusters in k-means clustering

FIGURE 6 :

 6 FIGURE 6 : Sum of squares method for finding the optimal size of clusters in hierarchical clustering

FIGURE 7 :

 7 FIGURE 7 : Comparing clustering methods' objective function (market-share=50%)

FIGURE 8 :

 8 FIGURE 8 : Comparing clustering methods' computation time (market-share=50%)

TABLE 1 :

 1 Example with 4 requests (configuration)

	Request Travel distance demand nshare EPT LDT
	1	19	1	3	8:00 8:45
	2	11	2	2	8:00 8:25
	3	24	2	1	8:15 9:00
	4	18	2	3	8:30 9:20

TABLE 2 :

 2 Example with 4 requests (solutions)

	Solution Total waiting time (s) Total travel distance (m) Number of cars
	A	28	3030	4
	B	51	2085	2
	C	53	1785	2
	D	29	2265	3

TABLE 3 :

 3 Clustering methods comparison

	Method Number of requests Objective function (normalized value) Computation time (s)
	Exact
	112	29.99	12966.00
	1092	1924.00	288000.00
	K-means clustering
	112	30.96	100.10
	1092	1994.18	1134.00
	4482	7761.03	5340.00
	11160	20773.70	20981.10
	Hierarchical clustering
	112	30.96	99.03
	1092	1994.59	1131.10
	4482	7791.90	5187.50
	11160	20905.87	19950.00

  1,092 requests. In both clustering methods, we try to have clusters with 50 requests to ensure potential trips to be shared inside the clusters for both test cases with 50 requests and 1092 requests. For 11,160 requests, the objective function is 0.64% lower for the k-means method, while the computation time is 5% more.In our Lyon6 + Villeurbanne test case, the maximum number of requests is 11,235. So we can use k-means clustering for this test case to have better solutions. However, for the Lyon network, we have more than 200,000 service requests. So the hierarchical method is a better solution for this network, since it is faster as it works directly with the similarity matrix, and it can provide high-quality solutions.

	a lower time.
	K-means clustering and hierarchical clustering methods increase the objective function by
	3.22% for 112 requests and 3.64% and 3.66% for 1,092 requests, compared to the optimal solu-
	tion. Both clustering methods can decrease the computation time from 216 minutes to less than 2
	minutes. This shows that both clustering methods are very effective in terms of reducing the com-
	putation time while keeping the quality of the solution acceptable. Then, by increasing the number
	of requests, the computation time for both methods exponentially increases (the major part of the
	k-means method computation time is dedicated to the multidimensional scaling method, which
	exponentially increases by increasing the size of the problem).
	K-means can give smaller objective function while hierarchical computation can result in
	TRB 2022 Annual Meeting	Original paper

TABLE 4 :

 4 Sum of Squares method for k-means clustering

	Number Size of clusters
	1	10
	2	20
	3	30
	4	40
	5	50

TABLE 5 :

 5 Sum of Squares method for hierarchical clustering

	Number Size of clusters
	1	75
	2	100
	3	125
	4	150
	5	175
	6	200
	7	225
	8	250
	9	275
	10	300

clustering methods on the shareability matrix to cluster the requests. Then, we described a method to solve the matching problem inside each cluster. To evaluate the method, we employed the proposed method on the network of Lyon city in France with half-million requests in the morning peak from 6 to 10 AM. The results showed that the proposed clustering method produces highquality solutions close to the optimal situation while reducing the computation time significantly.

The clustering method proposed in this study can be used for other similar systems, such as taxis and dial-a-ride services.

We have used earlier versions of this method to assess the impact of real-time ridesharing on road transportation networks (Alisoltani et al.,42). In our future works, we plan to assess its impact on large-scale multimodal networks, i.e. considering all available modes of transportation, and not only private cars using GeoTwin multi-agent simulation platform (Bauguion et al.,43).

AUTHOR CONTRIBUTION STATEMENT

All the authors have contributed to all aspects of this study, ranging from the conception and design of the methodology, analysis and interpretation of the results and discussion, and the manuscript preparation.

CONFLICT OF INTEREST

The authors declare no potential conflict of interests.