Hani Cheikh Sleiman
email: hani.cheikh.sleiman@outlook.com

Alessandro Tengattini

Matthieu Briffaut

Bruno Huet

Stefano Dal Pont

Hani Cheikh

Drying of mortar at ambient temperature studied using high resolution neutron tomography and numerical modeling

Introduction

The naturally occurring drying of cement-based materials has inevitable and significant effects on their durability. After casting, cement-based materials generally have a higher hygrometry than the surrounding environment. The water located in the porous network tends then to migrate towards the external environment to establish hydric equilibrium. This water loss induces drying shrinkage which gets restrained not only by the presence of aggregates, but also by the the core of the sample which remains saturated for a longer time. These differential deformations generate high stresses that can lead to cracking.

Additionally, the evolution of the relative humidity in the interconnected pores affects the hydraulic properties of the material, such as hydraulic conductivity and diffusivity [START_REF] Mualem | A new model for predicting the hydraulic conductivity of unsaturated porous media[END_REF] [START_REF] Perre | Simulation par volumes finis des transferts couplés en milieux poreux anisotropes: Séchage du bois à basse et à haute température[END_REF]), chemical properties (such as hydration degree [START_REF] Sciumè | A multiphysics model for concrete at early age applied to repairs problems[END_REF] and the kinetics of degrading chemical reactions [START_REF] Pan | Modeling of alkali-silica reaction in concrete: A review[END_REF]), thermal properties (such as thermal conductivity [START_REF] Baggio | Trasporto di calore e massa nel calcestruzzo cellulare[END_REF]), mechanical properties (such as creep [START_REF] Bazant | Concrete creep at variable humidity: constitutive law and mechanism[END_REF]), fire resistance [START_REF] Dauti | Analysis of moisture migration in concrete at high temperature through in-situ neutron tomography[END_REF] [START_REF] Tengattini | Quantification of evolving moisture profiles in concrete samplessubjected to temperature gradient by means of rapid neutrontomography: influence of boundary conditions, hygro-thermalloading history and spalling mitigation additives[END_REF], etc. It follows that characterising the hydric state and its evolution is pivotal for understanding its implications on the structural integrity of cement-based materials.

Several methods have been adopted in the literature to characterize the cement based material hydric state and properties. They can be broadly classified in three main categories: bulk-averaged or pointwise methods, full-field measurement techniques and indirect numerical modeling methods. Within the first group, the measure of the evolution of mass loss along time of a sample subject to drying is perhaps the most adopted method [START_REF] Baroghel-Bouny | Caractérisation microstructurale et hydrique des pâtes de ciment et des bétons ordinaires et à très hautes performances[END_REF] [START_REF] Baroghel-Bouny | Characterization and identification of equilibrium and transfer moisture properties for ordinary and high-performance cementitious materials[END_REF]. Despite its widespread adoption, perhaps related to the ease of its implementation, this measurement is a macroscopic average of the highly heterogeneous response of the materials and cannot provide insight about the spatial distribution of water content. To overcome this limitation, several authors have employed sensor-based measurement tools to assess relative humidity in a particular point in space [START_REF] Kim | Moisture diffusion of concrete considering self-desiccation at early ages[END_REF] [START_REF] Ma | Optimization of experiment methodology based on identification of parameters in concrete drying[END_REF]. However, these point-wise approaches are usually invasive and might show variable measurements depending on their proximity to different species (porosity, aggregates) in a heterogeneous porous medium such as concrete. Furthermore, they require the identification of the sorption/desorption isotherms of the material to relate the relative humidity to saturation [START_REF] Carette | Identifying the mechanisms of concrete drying: An experimental-numerical approach[END_REF].

In the last few decades, the development of full-field measurement techniques have allowed unprecedented insight into local processes [START_REF] Viggiani | Advanced experimental techniques in geomechanics[END_REF] [START_REF] Viggiani | Recent developments in laboratory testing of geomaterials with emphasis on imaging[END_REF]. When studying moisture transfer in porous media, multiple full-field techniques have been explored, such as Nuclear Magnetic Resonance NMR [START_REF] Faure | An MRI-SPI and NMR relaxation study of drying-hydration coupling effect on microstructure of cement-based materials at early age[END_REF] [17] [START_REF] Mcdonald | Sorption, anomalous water transport and dynamic porosity in cement paste: a spatially localised 1H NMR relaxation study and a proposed mechanism[END_REF], Radiation attenuation techniques (such as Gammadensimetry [START_REF] Villain | Gammadensimetry: A method to determine drying and carbonation profiles in concrete[END_REF]), Electrical techniques (Electrical Capacitance Tomography ECT [START_REF] Voss | Imaging moisture flows in cement-based materials using electrical capacitance tomography[END_REF]) and Electromagnetic Non-Destructive Evaluation ENDE techniques [START_REF] Loubser Du Plooy | Electromagnetic nondestructive evaluation techniques for the monitoring of water and chloride ingress into concrete: a comparative study[END_REF].

A particularly effective tool is Neutron Imaging, which lend itself as a very powerful non-destructive and non-invasive technique to study the moisture distribution in porous materials due to the high attenuation of neutrons by hydrogen [START_REF] Perfect | Neutron imaging of hydrogen-rich fluids in geomaterials and engineered porous media: A review[END_REF][23] [START_REF] Tengattini | Neutron imaging for geomechanics: A review, Geomechanics for Energy and the Environment[END_REF]. The 2D/3D attenuation fields resulting from the radiography/tomography can be correlated to water content after due calibration [START_REF] Stavropoulou | Dynamics of water absorption in callovo-oxfordian clayrock revealed with multimodal x-ray and neutron tomography[END_REF] [START_REF] Cheikh Sleiman | Simultaneous x-ray and neutron 4d tomographic study of drying-driven hydro-mechanical behavior of cement-based materials at moderate temperatures[END_REF].

Beside the intrinsic interest of understanding the evolving moisture profiles, the spatial distribution of moisture content is essential to rigorously calibrate thermo-hydric numerical models [27] [28]. This numerical calibration can, in turn, be used to characterize the hydraulic properties of the material which are anomalously characterized experimentally, such as the intrinsic permeability. In fact, the experimental measurement of this property shows a high sensitivity to the used fluid. The gas permeability can be up to three orders of magnitude higher than water permeability both because of pre-drying induced cracking and pore-wall physically adsorbed water removal. On the other hand, the water permeability can only be measured in saturated conditions, thus in a configuration where drying cracks effect on the drying kinetics are not taken into account. For these reasons, it is generally agreed that indirect numerical methods are more reliable than the aforementioned experimental ones in measuring the intrinsic permeability. However, due to the high number of model parameters and the mathematical flexibility of the constitutive laws, multiple solutions can be found, specifically when experimental data-sets related to local water content distribution are not measured [START_REF] Carette | Identifying the mechanisms of concrete drying: An experimental-numerical approach[END_REF].

In this contribution, an experimental/numerical approach is proposed to characterize the hydric state and deduce the hydraulic properties of a cement-based porous medium, taking advantage of the 3D nature of the acquired neutron tomographic data sets. In section 2, the experimental campaign is described, including the sample composition, their mass loss history and the neutron tomography campaign. A number of image processing operations and a novel methodology to estimate water content based on multi-phase porous media considerations are then presented in section 2.4 and 3. This in turn allows the characterisation of the saturation gradients at different instances in time. The resulting hydric gradients are then validated against the macroscopic mass loss history provided by gravimetric experimental data sets in section 3.3.2. These hydric gradients are then used together with the mass loss measurements as a complete data set to calibrate a common numerical model in section 4 by means of a minimisation algorithm detailed in section 4.3. Once calibrated, the model allowed an indirect assessment of the intrinsic permeability and the coefficient of hydric exchange.

Materials and methods

The applied procedure consists of scanning multiple mortar samples having different drying period in a thermo-hydrically controlled chamber using neutron tomography. Then, a series of image processing is applied to characterize the moisture profiles in these samples. The following subsections describe the samples preparation, the drying experimental apparatus, the neutron acquisition settings and the set of image processes that were used to transform grey scale fields to moisture content fields.

Specimens

The experimental campaign presented here involves 14 samples, all made out of EN standard mortar.

The water to cement ratio w/c is equal to 0.5 and the sand (standard sand CEN EN 196-1) to cement ratio s/c is equal to 3. These cylindrical samples were all cored out from a bigger bloc, which, after casting, was kept for two days in its mold in sealed conditions and was then placed in a bath of water to cure. A diamond coring bit with 12 mm inner diameter was then used to core out the samples. While one side of the samples was smooth due to its contact with the mold, the other side had to be further cut and smoothed with a rotating blade system, which was equipped with a water jet to avoid local thermal gradients. The resulting cylindrical samples had a height of 40 mm and were kept in water the beginning of the individual drying process to avoid any potential water exchanges with the environment. However, due to time constraints related to the beamline access, only 11 of 14 samples were chosen to be scanned based on their different drying periods.

Mass loss history

In order to induce different hydric states in the same material, the samples were placed at regular intervals in a drying environment. Each sample was first weighted and then rolled in three layers of adhesive aluminum tape on the longer side to avoid lateral moisture leakage and to induce a unidimensional moisture flow. The mass of aluminum tape was also measured for each sample to correct the measured mass loss. The samples were then put sequentially to dry in a thermo-hydrically controlled environment, with no ventilation, at In addition, two mortar samples were put to dry in an oven set to 105 • C to obtain a reference hydric state which is key for the quantification analysis as will be shown in the next section. Given the small size of the samples, the high temperature of the environment in which they were placed, and the duration of drying of 11 days, the mass loss of these two samples is here considered to correspond to the hydric equilibrium at 105 • C.

Neutron acquisition

Similarly to x-ray, neutron imaging requires an incident beam (I 0 ), emitted by a source (a nuclear reactor in this case) which is collimated through a pinhole. After the interaction with the sample the beam (I) is absorbed by a scintillator, which emits visible light proportionally to the incident beam intensity. This is focused and guided by optical elements and by a sCMOS camera. A series of these radiographies (2D attenuation fields) is acquired at different angular positions of the sample. These can then be reconstructed into a 3D attenuation field e.g., through a filtered backprojection method. The neutron tomographies in this study were acquired at the NeXT-Grenoble beamline (located at the Institut Laue-Langevin (ILL) in Grenoble, France and detailed in [START_REF] Tengattini | Next-grenoble, the neutron and x-ray tomograph in grenoble[END_REF]), which takes advantage of a uniquely high cold neutron flux. A schematic view of the setup is illustrated in Figure 2.

Before transporting the samples to the beamline, each sample was weighted and sealed to prevent any water exchange with the environment. Their weight was re-measured before and after each scan to ensure that no (or negligible) drying occurred during tomographic acquisition, which was found to be the case. An identical set of neutron acquisition parameters was selected for all the tested specimens. In absence of rapid dynamic processes occurring in the samples, these parameters were optimised to maximise spatial resolution.

In this case, 1088 projections were acquired in 2 hours and 54 minutes over a 360 • rotation, each radiography being the average of three acquisitions, each with 3.2 seconds exposure. To achieve this, an infinity corrected optical setup was employed, coupling a 50 mm Canon Lens with a 122 mm Heliflex lens, as detailed in [START_REF] Tengattini | Next-grenoble, the neutron and x-ray tomograph in grenoble[END_REF].

This setup yields a 16 µm pixel size. Given the small pixel size, a 20 µm thick Gd 2 O 2 S : T b 6 LiF scintillator was employed, A 30 mm pinhole D (corresponding to a neutron flux of 3x10 8 ncm 2 s at the sample, placed 10 meters away) was selected to achieve a true resolution of about 30 µm. Each neutron radiography was the average of 3 individual projections each with a 3.2 s exposure, resulting in a signal to noise SN R = 14.

Each set of radiographies was reconstructed into the corresponding 3D volume by means of the Feldkamp (FDK) back projection algorithm, as implemented in the commercial software X-act (from RX-Solutions).

An example of a 3D neutron attenuation field of the mortar sample is reported in Figure 3.

Segmentation, labelling and particle size distribution

Moisture transport happens essentially in the interconnected porosity of the cement given the comparatively low permeability of aggregates and sand particles. Therefore, in order to isolate the contribution of this phase, it is beneficial to identify and remove these particles from the neutron tomographies. Their identification (a process named segmentation in image processing) can in this case be achieved thanks to the lower neutron attenuation (i.e., the gray value of the voxels comprising them) of aggregates with respect to that of the cement paste. An histogram-based segmentation of the grains can then be achieved based on whether a voxel has a value above or below a suitably chosen threshold. In this work, to ensure repetability and avoid user dependency, this threshold was determined based on an algorithm named "IJ IsoData" [START_REF] Rider | Picture thresholding using an iterative selection method[END_REF] (available in ImageJ software [START_REF] Schindelin | Fiji: an open-source platform for biological-image analysis[END_REF]) that provides increasingly cleaner extractions of the object regions by means of an iterative process. This method allows the identification of the aggregates from the surrounding mortar (cement paste and sand fines below the imaging resolution) as illustrated in Figure 3. Once the particles are isolated it is straightforward to "remove" them from the subsequent elaborations. In order to quantify for the aggregates fines that were not captured by the neutron scans resolution, a granulometric analysis was done by applying the "watershed" morphological operation [START_REF] Beare | The watershed transform in itk-discussion and new developments[END_REF] on the segmented 130 image. This operation allows the separation of segmented aggregate phase into the individual grains comprising it as illustrated in Figure 3. The volume of each grain can then be measured and converted into a radius of a volume-equivalent sphere. This allows to calculate a particle size distribution as shown in Figure 4, where this cumulative probability density function is compared to the experimental sieving measurements.

Low attenuation

This comparison shows a good agreement above a 600 µm particle size.

Vertical gradients

A first indication of the evolution of the moisture content due to drying is reported in Figure 5, which shows the average grayvalue over the half-height of each sample (the half being completely within the field of view, given that the sample was larger than the field of view under for the chosen resolution). Each point From Figure 5, it can be observed that the reference samples with no mass loss (M15 and M16) show a slightly different initial water distribution. While in sample M15 the gray value attenuation profile suggests a uniform moisture content along its height, sample M16 shows an apparent moisture concentration near the top surface. This might be related to local variations in water content in between the mortar blocs from which the cylindrical samples were cored. Samples M15, M7 and M8 were in fact extracted from a different bloc different than the other samples. For this reason, in the upcoming analyses, the gray value of each sample is normalised to the corresponding reference sample cored from the same block. This measurement shows again the importance of the full-field nature of the proposed approach over the local, point-wise ones which could have led in this case to misleading conclusions.

Radial gradients

As aforementioned, the experimental campaign was conceived to impose a uni-directional drying front by sealing the lateral surface with multiple layers of aluminium tape. The three-dimensional nature of tomography allows the verification of this assumption through the study of the variation in attenuation values along the radial direction. To have a sufficient number of voxels for statistical representativity, the average over concentric (hollow) 20 slices-thick rings was compared. Additionally, the radial width of these sections was varied to maintain an almost-identical number of mortar voxels in each subvolume, as schematised in Figure 6a. The relatively small variation to voxel counts is due to the heterogeneity in aggregate content (which were removed from these statistics as detailed above). Voxel counters (red curves in Figure 6) are employed to verify that these small variations in voxel counts don't affect the gray value profiles across the rings. This operation was repeated at three different heights representing the top, middle and bottom of the sample, as reported in Figure 6 b, c and d, respectively. In each of these subfigures this operation is repeated for thee representative samples at three distinct drying stages. Figure 6 shows virtually no radial gradient between the core and the external surface of the sample, in all conditions and along the height. these results are in agreement with existing previous observations in [33][34], where adhesive aluminum was found to be effective at low temperature for periods shorter than 6 months. The samples were therefore considered to be perfectly sealed, justifying the hypothesis of unidirectional drying adopted in the following sections, which will only focus on vertical gradients averaged over the sample width.

(a) (b) (c) (d)

From attenuation values to saturation gradients 3.3.1. Saturation profiles

Once the aggregates are removed as detailed in section 2.4, the remaining material consists of partially saturated mortar (cement paste and sand fines below the imaging resolution). This can be further subdivided in three species: dry mortar "dm" (mixture of solid gel, unhydrated cement, chemically-bonded water and sand fines), capillary water "cw " and air "a".

The gray value (measured attenuation) of the mortar, considered as a multiphase mixture "mix " can be therefore expressed as the sum of these individual components:

µ mix = W cw µ cw + W a µ a + W dm µ dm , (1)
where W i and µ i are respectively the volumetric fraction and the intrinsic gray value of the specie i.

Assuming that the dry mortar undergoes a negligible mass variation (i.e., negligible effect of the dehydration reaction) and neglecting any contribution to the gray value induced by the modest volumetric variation (e.g., drying shrinkage), the variation of the compound attenuation between two saturation states can then be expressed as:

∆µ mix = ∆W cw µ cw + ∆W a µ a , (2)
where the increase in the volumetric fraction of air is equal to the reduction of water (∆W a = -∆W cw ) since the moisture leaving the porosity is replaced by air. The variation of capillary water volume fraction writes then:

∆W cw = ∆µ mix µ cw -µ a (3)
The gray value of air µ a can be directly measured in the tomographies and was found to average 7900 (gray value) with a standard deviation σ a of 1000. On the other hand, the gray value of water µ w was measured in multiple water-filled macro pores in the reference samples. It averaged 50500, with a standard deviation of 2000. It should be noted, nonetheless, that this value can slightly vary depending on the hydric state of the material, because of the scattering contribution to the attenuation. The attenuation law of water can in fact diverge from the beer-Lambert law, as generally the case for highly scattering elements in high concentration. For the specific energy spectrum of the used NeXT beamline, this discrepancy from the Beer Lambert was measured, for example, in [START_REF] Stavropoulou | Liquid water uptake in unconfined callovo oxfordian clay-rock studied with neutron and x-ray imaging[END_REF]. A possible approach to correct for this spurious fluctuation is then to account for a correction based on the observed correlation between equivalent water thickness "x" and water attenuation contribution, as measured on a reference, water-filled wedge sample. For thicknesses between 0 mm and 4 mm this can be fit with a second order polynomial µ w (x) = ax 2 + bx + c, where for this beamline a = -0.0052, b = -0.0023 and c = 0.3107. Assuming that this correction, measured in radiography, is also valid in tomography, the water gray value varies then non-linearly between 50500 and 57000 in between the reference, wet state and and ultimate, oven dry state.

The estimation of a saturation profile from capillary water volumetric fraction can be expressed as:

S l = S i l -∆S l = S i l - ∆W cw φ m (4)
Where S l is the saturation to liquid water, S i l is the initial, reference saturation, i.e., the early-age saturation after correcting for self-desiccation that reduces the fresh concrete relative humidity. φ m is the average porosity of the mortar phase after removal of the larger particles, as detailed above.

The volumetric fraction of fines ξ f below the resolution of the approach is estimated to be 40%, based on the granulometric analysis in section 2.4. thereby, the averaged porosity of the mortar after segmentation (cement and sand fines) can be estimated by interpolating between cement paste porosity and mortar porosity (with full granulometry), the interpolation coefficient being the fines volume fraction ξ f , it follows that:

φ m = φ cp -ξ f (φ cp -φ m ) (5)
where φ m is the mortar porosity, equal to 16%, and φ cp is the cement paste porosity estimated to 37% as reported in [START_REF] Rougelot | Etude multi-échelles des couplages hydriques, mécaniques et chimiques dans les matériaux cimentaires[END_REF] for the same types of material.

To finally deduce the saturation in Eq. 4, the only remaining parameter to determine is the initial saturation. To deduce this initial value, which can fluctuate depending on the exact sample history, the final, oven dried state (which is instead far more stable) was used as a reference. It should in fact be noted that variations in the choice of the initial saturation shift the ∆S l profiles to lower or higher absolute hydric states but do not affect the overall water loss along the tests. The initial saturation was then chosen to be 78% for all samples (despite the possible fluctuations between them), so as to ensure a saturation of less than 10% in the oven-dry states, as highlighted in Figure 7. This 10% limit was measured in [START_REF] Poyet | Experimental investigation of the effect of temperature on the first desorption isotherm of concrete[END_REF], [START_REF] Poyet | Temperature dependence of the sorption isotherms of cement-based materials: Heat of sorption and clausius-clapeyron formula[END_REF] in similar cementitious materials adopting a combined experimental and thermodynamics-based theoretical characterization of desorption isotherms up to moderate temperatures. The saturation profiles in Figure 7 show an overall gradient in between the upper exchange surface and the core of the sample up to 10%. This relatively weak saturation gradient is expected, given the lack of ventilation in the drying controlled environment. The additional local gradients along the samples' heights suggest local saturation variability in the porous networks of the samples.

Validation of mass loss measurements

The robustness of the grayvalue-saturation transition can be verified comparing it to the overall mass loss measured macroscopically over the whole sample via gravimetric analysis. While lacking the spatial information, in fact, this measurement is robust, and can be compared to the integral of the mass loss along the height of each sample. Specifically, mass loss can be calculated from each saturation profile S t l measured after a drying-period t and its difference from the reference initial saturation profile S 0 l as:

∆M t = ρ l . φ m l 0 (S 0 l -S t l ) ρ sat m , (6)
where ρ l is the water density, φ m is the averaged porosity of the mortar after removing the bigger particles as detailed above and ρ sat m is the density of mortar accounting for the full granulometry. The so-deduced mass loss for all the scanned samples, compared to the balance measurements is reported in Figure 8. The error bars for the imaging-based approach, show the effect of the aforementioned variability in the coefficient of attenuation of water and air. Despite the specificity in the composition of each sample, which can induce hydric variability as well as all the assumption previously stated, the gravimetry-based and tomographybased measurements are in good agreement.

Modeling of moisture transfer

To complete the characterisation of the processes at play, numerical models are of key importance to identify crucial hydraulic parameters, such as the intrinsic permeability and the coefficient of hydric exchange.

It is in fact possible to identify those parameters once the model is calibrated not only to the macroscale mass loss dataset but also the spatio-temporal gradients made available by the presented experimental campaign.

Several numerical models have been proposed in the literature to describe moisture transport in unsaturated porous media. They can be broadly classified into two main categories: complete and simplified models. While models belonging to the former establish a rigorous description of all the species of interest (liquid water, vapor and dry air) [39][40], models belonging to the latter propose a set of assumptions that simplify the physical description of the involved processes. This set of simplification generally remains valid when the gas pressure in the interconnected porosity remains equal to the atmospheric pressure, which means that all local pressure buildups dissipate quasi-instantaneously and the liquid permeation remains the dominant water transport mechanism. In this work, given the drying conditions imposed in the controlled environment (drying temperature of 20 • C and relative humidity of 35%), and given the material adopted (mortar having a water to cement ratio of 0.5), we assume the validity of a simplified model [41] [28].

Numerical model

The simplified model adopted derives from the mass conservation equation of the liquid water in the porous system, expressed as:

∂m l ∂t + ∇ • (m l v l-s ) = 0, ( 7
)
where m l is the mass of liquid water and v l-s is the velocity of the liquid water with respect to the solid skeleton of the material. The mass of water can be expressed in terms of density ρ l , porosity φ and liquid saturation rate S l as:

m l = ρ l S l φ (8)
whereas its velocity with respect to the solid skeleton can be computed using Darcy's law:

v l-s = -K k rl µ l ∇p l , (9)
where K, k rl and µ l are the intrinsic permeability of the system, the relative permeability of liquid water and its dynamic viscosity, respectively. Assuming the incompressibility of water, (ρ l = cte), as well as the existence of a differentiable function between saturation and capillary pressure, the conservation equation reads:

φ ∂S l ∂p c ∂p c ∂t -∇ • k rl (S l )K µ l ∇p c = 0 ( 10
)
This equation requires the definition of four quantities: porosity φ, a differentiable function of saturation, function of the capillary pressure S l (p c ), as well as intrinsic permeability K and relative permeability k rl .

Assuming that the hydration reaction of cement ended before the samples were exposed to drying, the microstructural properties, and in particular porosity and intrinsic permeability can be assumed as constant. The water retention curve (desorption isotherm) can be described by the Van Genuchten family of equations [START_REF] Van Genuchten | A closed-form equation for predicting the hydraulic conductivity of unsaturated soils[END_REF]:

S l (p c ) = - p c a 1 1-1 b + 1 -1 b , (11)
where the parameters a and b are fitting parameters. The Mualem [START_REF] Mualem | A new model for predicting the hydraulic conductivity of unsaturated porous media[END_REF] model is also used here to describe the relative permeability to liquid water as a function of water saturation:

K rl = S l 1 -1 -S b l 1/b 2 (12)
By analogy to the differential equation for heat transfer, it is possible to name φ ∂S l ∂pc the hydraulic capacity and Kk rl µ l the hydraulic conductivity. Equation ( 10) is here discretized and solved using the finite software Cast3m, and specifically its heat module accounting for the non-linear hydric conductivity term which depends on capillary pressure solution. The temporal discretization is done using the theta-method scheme.

Model inputs, initial and boundary conditions:

245

A quarter of the sample was described in axisymmetric conditions. A total of 182 eight-node elements were used for meshing the 2D surface, which was refined near the imposed boundary condition to ensure better convergence and higher results accuracy as shown in Figure 9. Two types of boundary conditions are generally employed in this kind of model to impose environmental conditions. The first is Dirichlet boundary condition, which consists of imposing the relative humidity of a surface (or an equivalent quantity, depending on the primary variable of the differential equation employed).

For example this could be used to impose the relative humidity within the environmental chamber directly at the exposed air-concrete interface. Nonetheless, this approach usually overestimates the interface exchanges, and neglects the presence of the so-called boundary-layer, whose moisture content evolves gradually depending on the distance to the material surface. The extent of this boundary-layer depends on the temperature of the chamber, the ventilation, as well as on the surface condition of the material. A second possible type of boundary conditions, are the Neumann-type boundary conditions, which can be generally expressed as:

-→ J bc = h bc .(p s c -p e c ). - → n , (13)
where p s c and p e c are the equivalent surfacic and external capillary pressures, respectively.

-→ J bc is the water flux and h bc is the convective exchange coefficient. This last parameter has a large variability (up to 4 orders of 250 0% 10% 20% 30% 40% 50% 60% 70% 80% 90% 100% magnitude), depending on the material and the experimental conditions [43] [44]. In the specific conditions of this experimental campaign, however, due to the low ventilation, it is expected to be small.

Relative humidity [%]

As mentioned in section 4.1, the model requires the definition of four quantities, namely porosity φ, a differentiable function of saturation with respect to capillary pressure S l (p c ), intrinsic permeability and relative permeability. All these properties were characterized experimentally in [START_REF] Rougelot | Etude multi-échelles des couplages hydriques, mécaniques et chimiques dans les matériaux cimentaires[END_REF], which has been used as a reference, with the exception of the intrinsic permeability (and more precisely, the permeability to water).

The initial relative humidity value is assumed equal to 80%, which corresponds to an initial saturation state of 78%, as deduced above. The correspondence between relative humidity and saturation is based on the desorption isotherm reported in [START_REF] Rougelot | Etude multi-échelles des couplages hydriques, mécaniques et chimiques dans les matériaux cimentaires[END_REF].

Model calibration and identification of hydraulic properties

Beside the aforementioned material properties, the model requires two additional fitting constants a and b, which influence the shape of the S l (p c ) and K rl (S l ) functions. Desorption isotherms for this class of materials were proposed in [START_REF] Rougelot | Etude multi-échelles des couplages hydriques, mécaniques et chimiques dans les matériaux cimentaires[END_REF], but, to the best of the the authors' knowledge, no direct experimental characterization of the relative permeability of the studied class of materials is proposed in the literature.

The parameter b, identified in the van Genuchten law, was therefore taken the same in the Mualem function for sake of consistency.

Once the isotherms are fitted against the dataset in [START_REF] Rougelot | Etude multi-échelles des couplages hydriques, mécaniques et chimiques dans les matériaux cimentaires[END_REF] and the consequent relative permeability function is identified, two parameters remain to be defined: the intrinsic permeability K and coefficient of hydric convection h bc . For this purpose, an optimisation procedure was implemented. It consists of minimising the error between the experimental and numerical response surfaces in the [K, h bc ] two-dimensional space. The goal is to find the values of K and h bc that minimize the total experimental/numerical error. This total error is composed of 3 normalized independent sub-errors defined as follow: the first sub-error computes the difference between the numerical mass loss output and the experimental one, averaged over the time steps.

The second sub-error computes the difference between the numerical saturation values and the experimental ones (previously shown in Figure 7), averaged over all the nodes on the geometric line [x 0 , x k ] (see Figure 9)

and then over all the time steps on which the measurements were done. Finally, the last sub-error computes the difference between numerical and experimental saturation "average slope", averaged over time. This "average slope" considers the saturation gradient as a line delimited by x 0 , where S x0 l = S max l , and x k , where S x k l = S min l . The system of equations of the numerical/experimental errors reads:

                                             tn ti=t0 ∆M t i exp -∆M t i num
Where ∆M ti exp and ∆M ti num are respectively the experimental and numerical mass loss at the time step t i , n is the total number of time steps, S l xj num and S l xj exp are respectively the numerical and experimental saturation at the node x j which belongs to the line [x 0 , x k ] going from the surface to the core of the sample. α, β, γ are the threshold factors under which the constraint is full field, specifically, α is the threshold factor 270 for the mass loss dataset, β and γ are respectively the threshold factor for the saturation value and the saturation slope data-set. The value of these threshold parameters are to be chosen by the user and depends on the tolerance he sets to achieve a satisfying fitting. In here, for instance, we assume that any mass loss error below 0.25% gives a satisfying fitting.

The key steps in the modelling approach and the experimental/numerical error minimization procedure The minimization procedure output is visualized in Figure 12. On the left, the evolving numerical and experimental mass loss surfaces in [K, h bc ] space are visualized at a specific time step, here, at 109 days. Each node of the numerical mass loss surface represents the mass loss output of the simulation having the corresponding [K, h bc ] as input parameters. As to be expected, this surface shows a minimum when h bc and K have there lowest value in the space, i.e., at K = 4 × 10 -21 and h bc = 0.1 and shows a maximum at K = 5 × 10 -20 and h bc = 3.5. On the other hand, the experimental mass loss response surface in this plot is represented by a constant horizontal surface in the [K, h bc ] space at each specific time step. This constant value corresponds to the interpolated mass loss value of the previously shown curve in Figure 1 at the corresponding numerical time steps. In the middle, the mass loss cumulative error response surface is shown together with the sets of [K, h bc ] that fulfills the first inequality in the equation 14. This set of points can be also represented by a polygon that joins the points on the boundary of the formed shape.

Finally, on the right, the graphical solution of each inequality in equation 14 is plotted. It was found that the intersection of these 3 polygons is a point that corresponds to a permeability K value equal to 10 -20 and to a coefficient of hydric exchange h bc equal to 0.4. The value of h bc was found to be one/two order(s) of magnitude smaller than reported in other studies (e.g., [START_REF] Soleilhet | Role of pore pressure on cracking and mechanical performance of concrete subjected to drying[END_REF]). This is likely caused by the lack of ventilation in the environmental control chamber. On the other hand, the value of K is marginally higher than the common estimation of this parameter. This discrepancy might be also attributed to the simplified nature of the adopted model which only accounts for the permeation of liquid water, neglecting the contribution of the diffusion of vapour. In fact, the low value of coefficient of hydric exchange slows the kinetics of the process, allowing a slow diffusion of vapor to take place in the porosity. The identified value of "K" is therefore the cumulative result of two distinct water transport processes: liquid water permeation and vapor diffusion. On the other hand, unlike other studies, the neutron tomographies acquired inhere allows for a direct measurement of saturation profile which is a more reliable dataset for the numerical fitting in comparison to other methods in the literature. Moreover, the suggested minimization procedure allows for a more objective fitting and for a more accurate identification for the hydraulic parameters than what can be found in different studies. Using the values of [K, h bc ] identified by the minimization procedure, the mass loss and saturation gradients numerical outputs are plotted against experimental data in Figure 13 and Figure 14, which shows a good agreement, within experimental error.

Conclusions

In this paper, an experimental-numerical approach was proposed to fully characterize the hydric state and properties of cement-based porous media. The novel aspect of this work lies in the adoption of neutron tomography to characterise the drying of cement-based materials at ambient temperatures with a major focus on measuring and validating water content profiles.

A set of 11 laterally-sealed cylindrical samples were dried in a thermo-hydrically controlled chamber for different time spans, while monitoring their weight loss along time. A high resolution neutron tomography was then acquired for each sample. The key phases were then differentiated in the images and their contribution to the overall neutron attenuation of the mixture was calculated. This allowed the estimation of the spatial distribution of water for different drying times with respect to a reference, wet sample. Based on an estimation of the initial saturation and the porosity of the material, it was then possible to deduce the saturation profiles to liquid water. These showed a gradient up to 10% between the exposed top surface and the core of the sample. The overall water loss of each sample deduced from the image-based saturation profiles and was also validated against the macroscopic mass loss derived from weight measurements.

After characterizing the spatio-temporal hydric state of the material, i.e. mass loss history and saturation gradients, a simplified finite element model was used to reproduce these experimental datasets and therefore identify key hydraulic parameters such as permeability K and coefficient of hydraulic conductivity h bc .

In addition, a minimization procedure was developed to ensure a more objective fitting procedure and to carefully assess the identification of those parameters. The numerical output were compared to the experimental ones and shows a fair agreement within the experimental error margin.

20 •

 20 C and 35 % relative humidity controlled with saturated MgCl 2 ,6H 2 O solution. The samples mass loss was characterised in time as shown in Figure 1.

Figure 1 :

 1 Figure 1: Drying of mortar samples: (a) Sample drying in a thermo-hydrically controlled environment (b) Total samples mass loss history

Figure 2 :

 2 Figure 2: A schematic view of the neutron instrumentation at the NeXT-Grenoble imaging beamline.

Figure 3 :

 3 Figure 3: Segmentation and labelling processing adopted in this work. To the left: rendering of a representative 3D reconstructed attenuation field. In the middle: rendering of the binary image after histogram-based segmentation . To the right: Labelled 3D volume isolating the individual grains, whose color corresponds to the assigned incremental index.

  granulometry of sand particles above 600 µm

Figure 4 :

 4 Figure 4: A comparison between sieving-based experimental cumulative granulometry and the corresponding one estimated through image analysis, starting from 600 µm to take into account the resolution of the neutron tomography scans

140

  in the plot corresponds to the average grey value mortar voxels in a 10 pixels-thick disk, along the height of the sample. The variation of gray value attenuation shown Figure5can be correlated to the variation of the hydric state along the height and drying stage as shown in the next section.

Figure 5 :

 5 Figure 5: Analysis of the attenuation profiles: (a) radiography showing the sample in the field of view (b) Grey value gradients of all the scanned samples.

Figure 6 :

 6 Figure 6: (a) Montage showing the partitioning of the sample in the radial direction so as to keep a similar number of pixels in each division. Evolution of grey value with respect to the radius for different saturation states and different positions in the sample: (b) Bottom (b) middle (c) Top

Figure 7 :

 7 Figure 7: Neutron-based measurement of saturation gradients of mortar samples at different hydric states. The measurement accounts for the variability of the fluids gray value attenuation coefficients.

  215

Figure 8 :

 8 Figure 8: Comparison between the balance-measured mass loss and the one deduced from image analysis.

Figure 9 :

 9 Figure 9: Scheme of the adopted finite element mesh and boundary conditions.

Figure 10 :

 10 Figure 10: The desorption isotherm and the relative permeability behavior laws required for the model closure.

Figure 11 :

 11 Figure 11: Flowchart detailing the modelling approach and the experimental/numerical error minimization procedure.

Figure 12 :

 12 Figure 12: Minimization procedure: (Left) numerical and experimental mass loss response surfaces in [K, h bc ] space at 109 days, (Middle) cumulative mass loss error surface together with the sets of [K, h bc ] fulfilling the mass loss inequality, (Right) Graphical solutions of all the inequalities in equation 14 which intersect at [K = 10 -20 , h bc = 0.4].

Figure 13 :

 13 Figure 13: Mass loss comparison of the: (red) balance measurements, (blue) image-based measurements and (cyan) FE estimation. The shaded area being to account for the variability of the acquired measurements.

Figure 14 :

 14 Figure 14: Saturation gradients: (shaded) image-based estimations accounting for the variability f the acquired measurement, (solid and dashed) FE estimations.

  for all 𝐾 𝑖 in 𝑲 for all ℎ 𝑏𝑐 𝑖 in 𝒉 𝒃𝒄 for all time steps in 𝒕 𝒏𝒖𝒎 • Solve the drying problem having 𝐾 𝑖 and ℎ 𝑏𝑐 𝑖 as input parameters • Compute the corresponding mass loss + export dataset • Compute the saturation profiles at 𝒕 𝒆𝒙𝒑 + export dataset Build the experimental mass loss surface response 𝑀𝐿 𝑒𝑥𝑝 𝐾, ℎ 𝑏𝑐 , 𝑡 = 𝑀𝐿 𝑒𝑥𝑝 (𝑡) • Build the experimental saturation volume response 𝑆 𝑙 𝑒𝑥𝑝 𝐾, ℎ 𝑏𝑐 , 𝑡, 𝑧 = 𝑆 𝑙 𝑒𝑥𝑝 𝑡, 𝑧 • Build the numerical mass loss response surface and saturation volume response • Compute the error between the experimental and numerical quantities • Plot the graphical solution of each inequality in the system (14) and find the intersection

	• Generate the mesh
	• Identify the parameters for desorption isotherm and relative permeability
	end
	end
	end
	Saturation profiles
	derived from neutron
	tomography at 𝒕 𝒆𝒙𝒑
	•

are illustrated in Figure11.

Data availability statement

The raw neutron tomography data-sets of the mortar samples will be available for ILL user club members in 3 to 5 years time at: https://data.ill.fr/proposal/getfile?file=201/d50tomo/exp_UGA-98/rawdata -the analyzed experiment are in folders 01 M1, 02 M7, 03 M11, 05 M14, 06 M15, 07 M6, 08 M16, 09 M8, 10 M4 and 12 M3.


Baptiste Gueuziec
email: baptiste.gueuziec@cea.fr

Jean-Pierre Gallois
email: jean-pierre.gallois@cea.fr

Frédéric Boulanger
email: frederic.boulanger@centralesupelec.fr

Qualitative Reasoning and Design Space Exploration

Keywords: Design Space Exploration, Qualitative reasoning, Constraint solving

The design of complex systems is a challenging task, combining optimization and testing techniques. Design space exploration allows the designer to optimize the parameters of the system, but is usually very costly and induces many inherent difficulties that require specific computation techniques to be solved. Qualitative reasoning was first introduced to describe system structures and causality, and developed to study the behavior of the system and discretize its state space into qualitative states. It is now used in diagnosis and verification to reduce computation time and precision while still preserving major properties of the behavior. This article presents how we think that qualitative reasoning can be applied to design space exploration in addition to state space discretization, i.e. how it may help in the choice of parameters for a system, by reducing the computational cost of this exploration.

INTRODUCTION

Qualitative reasoning is the area of computer science and engineering that focuses on the representation and reasoning on models with very little and imprecise information [START_REF] Forbus | Qualitative reasoning[END_REF]. Our previous works on this topic were applied to the study of hybrid systems and allowed us to study the system's state space based on its dynamics, constraints, and discrete controller. We obtained an entirely qualitative model of the behavior of a system, discretizing the state space into qualitative states that represent phases of the system's behavior. This discretization allows us to abstract both the state of the system and its behavior. The computation of the qualitative behavior generates path conditions requiring to be solved by symbolic solvers such as Sympy and Z3. In the instantiated systems, we obtained very satisfying results to study and anticipate the possible trajectories of the system. However, when we try to generalize our approach to the study of the design space, we meet the limits of these tools when computing the transitions. The design space of a system is the space of the possible values of its constant parameters, and its exploration requires more generality and abstraction than the exploration of the state space. Very few computational tools allow for the symbolic resolution of multivariate inequalities with variables of multiple orders. The difficulty comes from the fact that SMT solvers only find one solution that respects the given constraints, while the exploration of the design space would require finding all the values of the parameters for which there exists a solution in the variable space that satisfies the constraint. For a system with parameters P and variables X, exploring the state space amounts to seeking a valuation (A, x) of (P, X) satisfying a predicate R(A, x). When exploring the design space, we instead look for V P ⊂ D P such that ∀A ∈ V P , then ∃x ∈ D X verifying R(A, x) where D P and D X are the domains of the parameter P and the variable X. This logic problem cannot be solved generally by simply using classic SMT solvers.

Design Space Exploration (DSE) is a critical domain of the conception of complex systems. It consists in choosing, testing, and comparing different configurations of values for the parameters of the system before its construction to verify some property or optimize some utility function [START_REF] Aiguier | Complex Systems Design & Management[END_REF]. High-level DSE can be considered as multi-objective optimization [START_REF] Schafer | High-level synthesis design space exploration: Past, present, and future[END_REF] as we may want to optimize a set of conflicting constraints and utility functions. The solutions of the parameters are often searched in a finite space of authorized values [START_REF] Lattmann | Verification and design exploration through meta tool integration with openmodelica[END_REF], and the arising problem consists both in the explosion of the number of possible solutions when taking into account all the parameters and in the critical time needed to test the validity of all of these solutions. To this extent, an effective DSE requires three elements: a representation, an analysis, and an exploration method [START_REF] Kang | An approach for effective design space exploration[END_REF]. The representation of the system must be designed to optimize the research, the test, and the comparison of solutions. The analysis requires a metric, a utility function, or a set of logic predicates to evaluate and compare the different valuations of the parameters of the system. It also requires the choice of a verification method. For this verification task, various methodologies are employed, with different levels of reliability, depending on the complexity of the system and the number of solutions to test. These methods mainly consist of simulating the system to verify a property, measuring an output value, or using some verification process to check the validity of a specific property. Finally, the exploration method is necessary to choose the valuations to test and will strongly depend on the nature of the design space. For example, some contributions have studied the possibility of searching the design space using machine learning methods [START_REF] Blanchard | The uranie platform: an open-source software for optimisation, metamodelling and uncertainty analysis[END_REF][START_REF] Kim | Machine learning for design space exploration and optimization of manycore systems[END_REF][START_REF] Mudgal | Deep learning for entity matching: A design space exploration[END_REF], others use brute force search in discrete sets, and others use gradient descent methods in the case of continuous spaces [START_REF] Blanchard | The uranie platform: an open-source software for optimisation, metamodelling and uncertainty analysis[END_REF].

CONTEXT AND INITIAL PROBLEM

Hybrid Systems and Qualitative Modeling

We place ourselves in the context of hybrid cyberphysical systems, represented as:

Q the current mode of the system (discrete), X a set of n continuous variables, I a function mapping each mode to its invariant conditions, F a function mapping each mode to its flow conditions, T a set of discrete transitions, each described by their starting mode, their condition, their target mode, and their reset function P a set of m parameters Q takes values on D Q , the finite set of modes of the system, X is valuated on K n with K = R in the general case, and P takes values on D P a set of dimension m. D P can be either continuous or discrete. We make the realistic hypothesis that D P is at least bounded to fa-cilitate its exploration. We consider that the flow conditions are represented as ordinary differential equations (ODE), i.e., Ẋ = f (X,t) with t the time parameter. For example, an autonomous car running on a straight road on the axis x can be represented as:

D Q = {Stop, Accelerate, StableSpeed, Decelerate,
Reverse}, X = (x, ẋ) with x the coordinate of the car on the variation axis and ẋ its speed on this axis,

X = K 2 = R 2 , I = {Stop : X = 0 ; Accelerate : ẋ > 0 ; StableSpeed : ẋ > 0 ; Decelerate : ẋ > 0 ; Reverse : ẋ < 0}, F = {Stop : Ẋ = 0 ; Accelerate : Ẋ = X 0 0 1 0 + 0 a ; StableSpeed : Ẋ = X 0 0 1 0 + 0 0 ; Decelerate : Ẋ = X 0 0 1 0 + 0 b ; Reverse : Ẋ = X 0 0 1 0 + 0 0 } with a > 0, and b < 0 T = {(Stop, ẋ > ε, Accelerate, X → X), (Accelerate, ẋ = v max , StableSpeed, X → X), (StableSpeed, x > x max -γ, Decelerate, X → X), (Decelerate, ẋ < ε, Stop, X → (x, 0)), (Stop, ẋ < -ε, Reverse, X → X), (Reverse, ẋ > -ε, Stop, X → (x, 0))}
with ε > 0 a negligibility criterion under which a value is considered as null, γ a security distance from the objective point x max and v max the maximal speed allowed. P = {a, b, c, ε, γ, x max , v max } Qualitative modeling consists in the creation of a model adapted to support qualitative reasoning [START_REF] Sugeno | A fuzzy-logic-based approach to qualitative modeling[END_REF]. As this reasoning paradigm has the ambition to avoid numerical elements, the stake of qualitative modeling is to design models of Cyber-Phyical-Systems that integrate enough information from the system to perform efficient reasoning without using numerical knowledge or computation. We call qualitative model the model resulting from this design. A qualitative model must not consider variables as numerical valuations but as sets of values. For example, the autonomous car represented in subsection 2.1 can be qualitatively represented by imposing:

X = {-, 0, +}; F = {Stop : Ẋ = 0 ; Accelerate : Ẋ = ++ ; StableSpeed : Ẋ = +0 ; Decelerate : Ẋ = +-; Reverse : Ẋ = -0 }
In our previous works [START_REF] Gueuziec | Qualitative reasoning and cyber-physical systems: abstraction, modeling, and optimized simulation[END_REF], we studied the suitability of qualitative reasoning [START_REF] Kuipers | Qualitative simulation[END_REF][START_REF] Tiwari | Series of abstractions for hybrid automata[END_REF] to represent, study, and refine the state space of hybrid systems. We used polynomial solving to create a qualitative map of this state space and SMT solvers to determine all the possible transitions between the obtained qualitative states. This computation allowed us to generate a complete behavioral tree of the system between its modes and qualitative states. However, the limit of this approach is that it requires all the parameters of P to be known. Considering that one of the main causes of interest in qualitative reasoning is its high level of abstraction, its use in a context where all the values of P are not entirely defined seems appropriate. However, since the computation of the transitions in the state space of the system requires solving symbolic inequalities with a SMT solver (we use Z3 in our model), the presence of symbolic constants in the equations of the system creates significant problems. These symbolic constants should not be treated on the same logic level, and using a simple SMT solver as previously done is not sufficient anymore because we are now searching for a subset V P ⊂ D P of parameter values such that for any set p of parameters in V P , there exists x ∈ K n that satisfies a predicate R(p, x). For example, if we consider a Van der Pol oscillator system whose dynamics is given by the equation system :

ẋ = 10(y + x -x 3 3 ) ẏ = p -x -3y 4 (1)
with p a positive parameter, we can look for a value of p such that the system is convergent. To this extent, the expression of the predicate R(p, X) would be: ẋ = ẏ = 0, i.e. 10(y + x -x 3 3 ) = px -3y 4 = 0. The DSE implies to look for a value of p > 0 such that ∃(x, y) such that 10(y + x -x 3 3 ) = px -3y 4 = 0. SMT solvers require well-designed predicates on the qualitative states and frontiers and allow us to check composed properties containing conjunctions and disjunctions, meaning that it is either possible to operate on the state space by giving the satisfiability for a valuation (complete or partial) of P, or on the design space to evaluate the pertinence of a choice of P. Our motivation in this work was to apply qualitative reasoning and analysis to the simulation of cyberphysical systems. Therefore, our motivation is to develop our methodology to be able to use qualitative reasoning for DSE to choose or optimize the values of the parameters P to satisfy some constraints.

Design Space Exploration

Design Space Exploration aims to find a design satisfying a set of constraints Const or optimizing some utility function f . Depending on the situation, the formulation will have the form of a satisfiability problem or an optimization problem. Since we want to explore the possibility of using qualitative reasoning, we mainly focused on the satisfiability case. Therefore, the problem takes the form:

max(V p = (V 1 ,V 2 , ...,V m ) ⊂ D P ) s.t. Const(V p ) = True
where Const(V P ) represents the complete set of constraints that the parameters must satisfy. As qualitative reasoning is based on the loss of numerical values, optimization problems are not adapted for this modeling paradigm.

The sought maximum of the subspace V P is a maximum regarding inclusion. Getting the maximum valid set regarding the design constraints is essential to compare the allowed values on other criteria afterward. This comparison is at the base of the design optimization [START_REF] Fuchs | Discrete search in design optimization[END_REF] and requires further analysis steps. They will not be treated in this article.

The contribution presented here aims to study the possibility of linking the DSE problem with qualitative reasoning, as already mentioned in [START_REF] Lattmann | Verification and design exploration through meta tool integration with openmodelica[END_REF]. Our goal is to open the capabilities of qualitative reasoning to broader perspectives and to present a hint at its possible use for more complex applications than state space partitioning. Moreover, since DSE as presented in [START_REF] Kang | An approach for effective design space exploration[END_REF] is mainly limited to structural considerations, the addition of qualitative reasoning to the DSE tools may increase the field of resolution with dynamic and behavioral considerations.

QUALITATIVE REASONING

As explained in the introduction, the first significant element required to process the exploration of the design space is an adapted representation of the model. Using a discretization process, we get a qualitative representation of the system corresponding to QM = ⟨Q, X, B, F, I, T, P⟩ with B a mapping from each mode q in Q to a set of qualitative borders defining the qualitative state they separate.

Definition: We call qualitative state of a system the pair (m, qs) associated with a list of polynomials P with m ∈ D Q the current mode and qs ∈ {+, 0, -} |P| a vector such that ∀i ∈ 1, |P| , qs[i] =if P[i](X) < 0, 0 if P[i](X) = 0, and + otherwise.

The qualitative states are then represented as an array of {-, 0, +}, showing the position of the corresponding state with respect to each frontier in B. This representation can be applied to fully designed systems or more abstract ones with non-valuated parameters of P. The second case implies the impossibility of computing the transitions between qualitative states, as this knowledge requires the computation of the sign of the Lie derivative [START_REF] Tiwari | Series of abstractions for hybrid automata[END_REF] to deduce the allowed transition directions.

The Lie derivative [START_REF] Yano | The theory of Lie derivatives and its applications[END_REF] consists in a vector field differentiation method enabling the directional derivation of a function. For a polynomial p whose zeros define borders of the system, the Lie derivative corresponding to p can be written as

L X (p) = ∑ X i ∈X ∂p ∂X i ∂X i
∂t where the X i are the components of X and t is the time parameter. Determining the sign of this expression requires the ability to fully evaluate p and the derivatives of all the variables X i .

This computation allows the creation of a qualitative model, where the states and evolution are not represented with numerical values but with variable signs (+, 0, -) and variation directions (corresponding to the sign of the associated derivative).

As the DSE requires optimization criteria or necessary constraints, it is possible to integrate these elements in the construction of the frontiers during the abstraction process. If R(P, X) ∈ K n+m is a predicate on X and P that should be satisfied by the system, then the equations composing R(P, X) can be considered individually to define new borders. We anticipate that qualitative reasoning may improve the DSE efficiency and applicability to complex systems as qualitative reasoning applied on early stages of the design of a system can link the found solutions to non valuated variables by applying purely symbolic solving to estimate adapted values of P to different use cases of a same system depending on the anticipated variations in the state space evolution between different case studies. Moreover, the abstraction of the state space may allow to propose local solutions to adapt the design of the required qualitative behaviors.

PARAMETERS EVALUATION

The second important aspect that is necessary for DSE is an analysis method. When a value or a subset of values is available for P, it is essential to test its pertinence using a solver adapted for the equation structures. Here, as we made the hypothesis of systems defined with polynomial equations and predicates, we use a symbolic polynomial solver. Many behavioral properties of a dynamic system can be conveniently represented using temporal logic [START_REF] Han | Linear temporal logic for hybrid dynamical systems: Characterizations and sufficient conditions[END_REF]. However, predicate satisfiability for temporal logic is more complex and costly than first-order logic (FOL) evaluation. Moreover, qualitative reasoning gives us an adapted support to reason on firstorder logic predicates. We will limit ourselves to the properties expressible as FOL predicates to apply the possibilities of qualitative reasoning. This choice implies representing the behavioral constraint we want to express as predicates.

For example, let us suppose we want a deterministic behavior for the system modeled through its discretization in qualitative states. In that case, we may wish to have a direct transition from the qualitative state (m 1 , s 1 ) to (m 1 , s 2 ). In temporal logic and considering a temporal structure T, <, we could write the predicate corresponding to a necessary transition q s = (m 1 , s 1 ) U q s = (m 1 , s 2 ) with s 1 and s 2 qualitative states of the system in mode m 1 , i.e., abstractions of values of X following an abstraction function adapted to the system. In FOL, we could write it using the Lie derivative such that

(m 1 , s 1 ) =⇒ ∀b ∈ B(m 1 , s 1 ) \ B(m 1 , s 2 ), L X (p b ) * p b (s 1 ) > 0 ∧ for b f ∈ B(m 1 , s 1 ) ∩ B(m 1 , s 2 ), L X (p b f ) * p b f (s 1 ) < 0 where B(m 1 , s i )
is the set of borders defining the qualitative state s i and p b the polynomial function whose zero corresponds to b. This predicate corresponds to the evaluation of the sign of each L X (p b ) and its comparison to the evolution of p b after crossing the corresponding qualitative frontier. A qualitative description of the system gives us a convenient toolbox to express the dynamic logic of the system as first-order logic predicates.

In our case, we keep using Z3 as it perfectly handles property verification on parameter values. We use it as a Python library and create a solver s to which we add the constraints on the system's design, as will be shown later in the algorithms we use for design space exploration. The command s.check() returns the satisfiability of the created problem under the given constraints. Compared to constraint solving on the state space, the main difficulty here is that there must be many calls to the solver to fix a value or a set of values of P that verify the constraints.

Finite Design Space

In the literature, DSE methods are mainly adapted to optimize the search for solutions in finite design spaces. Many practices have emerged, including genetic algorithms [START_REF] Palesi | Multi-objective design space exploration using genetic algorithms[END_REF], machine learning methods [START_REF] Motamedi | Design space exploration of fpga-based deep convolutional neural networks[END_REF], and other discrete optimization methods. We chose to exploit the permissiveness of the SMT solver Z3 to use branch-and-bound to specify the different parameters of P sequentially without dealing with every possibility. In branch-and-bound, we compute for each parameter p i of P the satisfiability of the constraints set for every proposed value v i, j associated with every satisfiable combination of v k,l for k < i. Here, v i, j corresponds to the j th authorised value of the i th element of P. By cutting the branch of the non-satisfiable partial combinations of P 1-i , we save the computation time that could be lost to test trivially impossible solutions. This algorithm is detailed in algorithm 1, where possibleCon f ig retains the authorized partial evaluation allowing the constraints Constr to be satisfiable. The function toConstr transforms the instantiation of the current parameter p i to a constraint expressed as an equality. This new constraint is noted newConstr, and its conjunction with Constr is noted compConstr. D P [p i ] contains all the authorized values v i, j of the component p i of P. The operation pre + v corresponds to the association of the authorized partial configuration pre of the p k , k < i with the considered value v of p i .

Infinite Design Space

Although the exploration of discrete design space is often practiced on finite discrete sets, developing the same techniques for continuous or infinite research spaces can be crucial. We develop an approach based on dichotomy, i.e., a successive partition of the search space in two distinct parts. For a single parameter p, we define a minimum width of the remaining definition set, and we subdivide it until we find that there is no solution in the set or that the current width of the set is smaller than the minimum width. The algorithm is shown in algorithm 2, where DichotomialCut(S p ) proceeds to a split of the set S p . In the general case, if S p is an interval on R, noted [a, b], the separation is completed by computing the center of the interval c = a+b 2 . Otherwise, if S p is composed of many convex non-contiguous sets, the separation is done by splitting S p between two noncontiguous components. minSize corresponds to the criteria on max(∥a -b∥) a,b∈set to stop the splitting process. As in algorithm 1, Constr represents the set Data: P, D P , Constr Result: Possible values of the parameters in P possibleCon f ig = [[]] * length(P) for i ∈ (0, length(P) -1) do

p = P[i] for pre ∈ possibleCon f ig do for v ∈ D P [p] do newConstrs = toConstr(pre + v) compConstr = newConstr ∧Constr s = SMT solver s.addConstraints(compConstr) an = s.Check() if an == SAT then possibleCon f ig[i + 1].append( pre + v) ) end if isEmpty(possibleCon f ig[i + 1]) ∧ i < length(P) -1 then Stop end end end end return possibleConfig Algorithm 1: Branch-and-Bound.
of constraints on the design parameters. allowedSets contains all the explored sets and Boolean values corresponding to the satisfiability of the design constraint on each of these sets.

To deal with multi-dimensional P, we can combine algorithm 2 and algorithm 1 to use dichotomy on every p i ∈ P sequentially and to divide the different combinations with a disjunction of cases.

Interestingly, we process once again to a discretization of a continuous space just as in the state space partitioning. The obtained partition of D P is based on the respect of the design constraints. Using an optimization function could also give more tools to refine such a partition based on its magnitude and variation. Moreover, the computation of its derivative would allow us to compute its critical points and, therefore, its maximums.

Non-Uniform Probabilistic Distribution

In the main part of the cases, D P is a subset of R and often takes the form of an interval. However, in some situations, it may appear that the design space does not take this form. For example, if prior knowledge is available regarding the probability of the different possible values, the set can no longer be represented using a simple interval. Some works, such as [START_REF] Blanchard | The uranie platform: an open-source software for optimisation, metamodelling and uncertainty analysis[END_REF], proposed solutions using probability density functions to represent this irregular distribution. This can be seen as a generalization of the non-probabilistic case considering that ∀(a, b) ∈ R 2 such that a < b, the interval [a, b] can be represented by a probability distribution following the uniform probability density law U(a, b) on the same interval.

One can also represent discrete sets using discrete probability laws. Therefore, it is possible to model many design spaces using probability functions. Using the expression of the probability density function for each p ∈ P, it is possible to pick values of P according to these density functions and combine this biased selection with the previously presented methods. However, this way of exploring the design space better fits the classical numerical simulation verification as it also applies gradient descent to orient the choice of new values for the next simulation. It constrains the differentiability of the utility function according to the researched element. This also better corresponds to the selection of the variable's initial values than to parameters valuation.

EXPERIMENTATION

At this stage of our research, we have not been able to apply the presented algorithm to concrete cases, and we did not set up a test benchmark either. We intend to implement this in the following weeks to put our propositions to the test. Specifically, we plan to test these algorithms on both simple systems such as a Brusselator system, with positive parameters a and b and state variables x and y, whose dynamics is given by Equation 2, and on more complex ones such as the Lorenz system, with positive parameters σ, ρ and β and state variables x, y, and z, whose dynamics are given by Equation 1.

ẋ = 1 -(b + 1)x + ax 2 y ẏ = bx -ax 2 y (2)    ẋ = σ(y -x) ẏ = x(ρ -z) -y ż = xy -βz (3)
Both systems exhibit different behaviors depending on the value of their parameters. The convergence or periodicity of such systems being an essential element in the trace of a system, we will try to express the sets of validity of the convergence of each one for finite and continuous sets of allowed values. If our method allows to obtain sufficiently precise results corresponding to theoretical calculations, it will be possible to consider further analysis and computation to improve the efficiency and the generality of the process. However, to consider the result interesting, it will be necessary to measure the computation time of our tests to ensure that there is either a gain in complexity, permissiveness, or generalization.

PERSPECTIVES AND FUTURE WORKS

The next step to perfectly combine our different subjects would be to express the properties of the state space and the solution of the various constraints on the variable X depending on the chosen value of the parameters P. The validity of a predicate would be expressed as X = f (P) for the subset of D P that allows these predicates to have a solution. To achieve this goal, we searched many different methods and looked at other tools, such as modal solvers or Wolfram Alpha, which give encouraging solutions in that direction. However, it does not generalize well for systems with too many variables. Other studied tools allowed this type of resolution but only for low-degree polynomial functions. There is no possible generalization for now.

Some works have studied the different possibilities offered by qualitative reasoning. [START_REF] Medimegh | Qualitative simulation of hybrid systems with an application to sysml models[END_REF] focused on the description and the prediction of the behavior of systems at a high level of abstraction, while [START_REF] Zaatiti | Abstractions refinement for hybrid systems diagnosability analysis[END_REF] developed the diagnosis aspect of qualitative reasoning using more numerical methods. [START_REF] Gueuziec | Qualitative reasoning and cyber-physical systems: abstraction, modeling, and optimized simulation[END_REF] proposed a technique to automatize system abstraction and the creation of hybrid automata. All these projects focused on the state space study and did not delve into state space partitioning but showed the advantages and the limits of the different approaches. Many works proposed very different approaches regarding design space exploration, from the classic search in a finite design space [START_REF] Lattmann | Verification and design exploration through meta tool integration with openmodelica[END_REF] to works presenting deep-learning-based strategies [START_REF] Motamedi | Design space exploration of fpga-based deep convolutional neural networks[END_REF]. Methods also exist to deal with continuous sets using probability density functions [START_REF] Blanchard | The uranie platform: an open-source software for optimisation, metamodelling and uncertainty analysis[END_REF] to represent the design space based on previous knowledge and more qualitative constraints regarding the expected value of parameters. However, this approach is more appropriate to deal with the choice of the initial value of variables and is, therefore, at the edge between state and design space exploration. Also, many works regarding the solving of temporal logic constraints and its use for optimization have been led [START_REF] Wolff | Optimization-based trajectory generation with linear temporal logic specifications[END_REF] and may allow significant progress in the resolution of temporal and modal predicates. DSE, in its aspect of design optimization, has been more significantly treated in [START_REF] Fuchs | Discrete search in design optimization[END_REF][START_REF] Wang | Review of metamodeling techniques in support of engineering design optimization[END_REF].

CONCLUSION

Theoretically, qualitative reasoning should have an interesting relation with DSE, and we can take advantage of it to optimize performances and computation time. Qualitative abstraction allows a permissive representation of models and a satisfying logical expression of essential system properties. Using an SMT solver on these first-order logic predicates is a relevant method to test the constraints' validity for a specific valuation of the parameter list P.

  Data: p, S p , Constr, minSize Result: Domain of possible values for p allowedSets = [] S b , S u = DichotomialCut(S p ) for set ∈ {S b , S u } do

	newConstr = toConstr(p ∈ set)
	ConstrTot = newConstr ∧Constr
	s = SMT solver
	s.addConstraint(ConstrTot)
	an = s.Check()
	if an == Unsat then
	allowedSets.append((set, False))
	else
	if max(abs(a -b)) a,b∈set < minSize
	then
	allowedSets.append((set, True))
	else
	allowedSets.concatenate(
	Dichotomial search(
	p, set,Constr, minSize
	)
	)
	end
	end
	end
	return allowedSets
	Algorithm 2: Dichotomial search

DESIGN SPACE EXPLORATIONThe last element of an efficient DSE is an optimized research method to choose the values or the sets of values of P. As already mentioned, we assume that the design space is bounded. To develop a research method, we have to separate two prominent cases depending on the nature of the Design space.