domain: phys sde

text:
Adrien Toesca

Damien David
email: damien.david@insa-lyon.fr

Kévyn Johannes

Michel Lussault

Generation of weather data for the assessment of building performances under future heatwave conditions

Keywords: Heatwaves, global warming, building simulations, weather data

Heatwave weather files are needed to assess building performance under future heatwaves. This paper presents a methodology for producing a minimum set of heatwave weather files, which should be representative of the diversity of heatwaves expected in a location of interest. It is a four-step methodology. Weather projections are first collected from the CORDEX project database. Then, heatwaves are identified in the weather data. A list of independent and significant indicators is constructed to

Introduction

During the next decades, it is highly probable that a majority of the buildings on earth will experience several heatwaves. Additionally, the heatwaves to which each building may be exposed are likely to have different characteristics. Some of them will last longer, some of them will have higher temperatures, etc. The wide variety of heatwaves to which a building could potentially be exposed can be attributed to the stochastic nature of extreme meteorological phenomena. It is also due to the uncertainty in the evolution of the climate, which may be attributed to both the uncertainty in the evolution of greenhouse gas emissions and the limitations of the models used to predict future climate evolution [START_REF] Machard | A Methodology for Assembling Future Weather Files Including Heatwaves for Building Thermal Simulations from the European Coordinated Regional Downscaling Experiment (EURO-CORDEX) Climate Data[END_REF]. Considering the diversity of the heatwaves that a building may have to face during the next decades, a question arises: is it relevant to assess the performance of the building during future heatwaves by simulating the building thermal behaviour under the conditions of a unique heatwave? To the authors' point of view, most of the time it is not relevant. Instead, a set of heatwave weather files would be needed to assess the building performance, with the heatwaves being representative of the diversity of the future probable heatwaves at the building location. Some authors have developed and tested different methods to constitute such a set of heatwave weather files. [START_REF] Ramakrishnan | Thermal performance of buildings integrated with phase change materials to reduce heat stress risks during extreme heatwave events[END_REF] studied the thermal performance of buildings equipped with phase change material panels under heatwave conditions in Melbourne. They used the weather data of the 2009 heatwave as a base file. To consider the impacts of climate change, they constructed heatwave weather data for 2030 and 2050 by increasing the air temperature of the 2009 heatwave by 1 °C and 1.5 °C, respectively. They did not modify the other weather variables. This method does not produce heatwave datasets that are representative of the diversity of future heatwaves because the weather data are too similar. [START_REF] Lassandro | Multi-criteria and multiscale assessment of building envelope response-ability to rising heat waves[END_REF] studied resilient retrofitting strategies of buildings under future heatwave conditions. They used three heatwave weather files. Those weather files were produced by selecting the hottest summer week of the representative International Weather for Energy Calculation (IWEC) weather data of three Mediterranean cities located in different climate zones. The method relies on the assumption that the time variability of future heatwaves could be similar to the spatial variability of current heatwaves.

[4] developed a method for defining reference summer weather year (RSWY) weather files from historical weather data over a period of 31 years. The method defines a heatwave as an event during which the t-Standard Effective Temperature (t-SET) comfort index, computed with weather data, remains over an arbitrary threshold for several consecutive days. The authors calculate three indicators for each heatwave: the duration, the severity (surface between the t-SET time evolution over the threshold) and the intensity (ratio between severity and duration). The RSWYs are the second most extreme years (over the 31 historical years) in terms of heatwave duration, intensity or severity. There can therefore be 3 RSWYs for a given period and location. It is also possible that the same year is the second most extreme on all 3 indicators. To the present authors' knowledge, the method has only been applied to historical weather data. [START_REF] Yaqubi | Creation and application of future typical weather files in the evaluation of indoor overheating in freefloating buildings[END_REF] built Future Typical Weather Year (FTWY) weather files from weather projection data for the location of Nantes, France. The authors compared the FTWYs obtained with weather data produced by three different climate projection models to the weather data of 2003, during which a deadly heatwave occurred in France. They observed that the heatwaves contained in their future typical weather data years were less extreme than the 2003 heatwave. This may be true for the heatwaves contained in typical weather files, but this may not be true for more extreme future heatwaves. In the authors' view, extreme weather files should only be compared with extreme files and not with typical weather files.

None of the methodologies mentioned here can be used to produce a representative set of heatwave weather files for building assessment. The methods for producing representative meteorological years do not focus enough on the heatwave characteristics. Most of the methodologies are only based on an analysis of the air temperature evolution, but a building's thermal behaviour strongly depends on other weather variables, such as solar radiation, wind speed, or sky temperature. Finally, none of the methodologies are able to assess the representativeness of the set of heatwave weather files regarding the possible diversity of future heatwave characteristics.

The present paper describes a methodology that intends to construct a minimum representative set of heatwave weather files for building assessment. It is a four-step methodology. An overview of the methodology is proposed in the first section. A detailed description of each step is provided in the following sections. The description of the methodology is illustrated by an application for the location of the Lyon-Saint Exupéry Airport.

Methodology overview

Figure 1 summarizes the overall methodology. The first step consists of collecting weather projections for the locality of interest. Here, the term "weather projection" refers to a long-term weather prediction. The second step consists of identifying heatwaves in the weather projections. This data produces the extended heatwave dataset. The third and fourth steps intend to sample the extended heatwave dataset. The third step consists of building a heatwave characterization space in which each heatwave is located according to its characteristics. The fourth step consists of selecting heatwaves in the characterization space. For the purpose of this paper, the four-step methodology has been applied to the case study of the Lyon_Saint Exupéry Airport, France. Table 1 provides general information about the case study. The Lyon-Saint Exupery Airport is located in a rural area, 30 km away from the city of Lyon, France. The method is applied in rural areas and therefore does not consider urban heat island effects. Thus, the results of the method will not depend on assumptions about the evolution of cities, which will allow the set of heatwaves to be used later as a base to simulate the urban heat island effects during future heatwaves for different urbanisation. The future heatwave weather data are intended to be used for the evaluation of the absolute building performances during future heatwaves, as opposed to a comparative evaluation with a reference building behaviour evaluated during a reference historical period. In that context, there is general agreement on the fact that bias-corrected weather data should be used [START_REF] Yaqubi | Creation and application of future typical weather files in the evaluation of indoor overheating in freefloating buildings[END_REF]. In the CORDEX database, bias-corrected data are available only for 3 variables: the air temperature, the wind speed and the global solar irradiation.

Only the weather projections for which bias-corrected data are available for those 3

variables are considered suitable.

The first step of the methodology consists of identifying suitable weather projections according to the three criteria presented above, downloading the related weather data, computing the missing weather variables, and interpolating the data linearly to obtain the full list of data that is required in building performance simulation weather files.

Application

Table 3 shows the result of the weather projection selection for the location of Lyon-St Exupery. The Lyon St-Exupery site is located in the EURO-CORDEX subdomain of the CORDEX project. The selection of weather projections (Table 3) is relevant for any location in this subdomain.

Fourteen weather projections were selected. They correspond to 2 RCP scenarios, RCP4.

CNRM-CERFACS-CNRM-CM5 RCA4 r1i1p1 X X ICHEC-EC-EARTH RAMCO22E r1i1p1 X X IPSL-IPSL-CM5A-MR RCA4 r1i1p1 X X MOHC-HadGEM2-ES RAMCO22E r1i1p1 X X MOHC-HadGEM2-ES RCA4 r1i1p1 X X MPI-M-MPI-ESM RCA4 r1i1p1 X X ICHEC-EC-EARTH RCA4 r12i1p1 X X 3.2.
Step 2: Identifying heatwaves

Objectives and method

The second step of the methodology consists of identifying heatwaves in the weather projections. There is no universal definition of a heatwave, but the literature frequently refers to four heatwave identification methods: the "Spic, Sdeb, Sint" method from [START_REF] Ouzeau | Heat waves analysis over France in present and future climate: Application of a new method on the EURO-CORDEX ensemble[END_REF],

the "IBM" method from [START_REF]France's heat health watch warning system[END_REF], the Excess Heat Factor ("EHF") method from [START_REF] Nairn | The Excess Heat Factor: A Metric for Heatwave Intensity and Its Use in Classifying Heatwave Severity[END_REF] and the HeatWave Magnitude Index ("HWMI") method from [START_REF] Russo | Magnitude of extreme heat waves in present climate and their projection in a warming world[END_REF].

To avoid inconsistencies in the heatwave dataset, a choice must be made between the four identification methods. It is not possible to base that choice on the description of the four identification mechanisms. Therefore, it was decided to perform heatwave identification with the four methods and to identify the most appropriate method from a comparison of the resulting heatwave datasets.

Application

Table 4 contains a statistical description of the four heatwave datasets for the location of Lyon-St Exupery. For each heatwave dataset, the first column contains the average number of identified heatwaves per year, the second column contains the average duration of the heatwaves, the third column contains the average value of the average temperatures during the heatwaves, and the fourth column contains the absolute number of identified heatwaves.

The "EHF" and "HWMI" methods identified a considerably larger number of heatwaves than the first two methods. Their identification criteria are less restrictive, with lower threshold values. As a result, the average temperature during a considerable part of the heatwaves is mild. However, health and discomfort issues might be expected during more extreme events. For this reason, those two methods were not chosen.

The distinction between the "Spic, Sdeb, Sint" method and the "IBM" method is more subtle. The two methods identified roughly the same number of heatwaves, and a deeper analysis of the datasets reveals that 90% of the heatwaves identified by the two methods correspond to the same extreme weather event: the same period of the same weather projection. The major difference between the 2 methods is the criterion used to locate the beginning and end of the heatwaves. This criterion is less restrictive for the "Spic, Sdeb, Sint" method in which the heatwaves start earlier and finish later. As a consequence, the heatwaves last for longer (8 days on average), with a lower average air temperature (2 °C on average). The "IBM" method only identifies the heart of heatwaves. Although those observations reveal some specificities of the two heatwave identification methods, they did not motivate the choice between the two methods.

Finally, the "IBM" method was chosen. The main reason for that choice is that the method is used by France's Heat Health Watch Warning System to trigger heatwave alerts [START_REF]France's heat health watch warning system[END_REF].

The second step of the methodology, applied to the Lyon St Exupery, resulted in the construction of a dataset containing 2229 heatwaves (highlighted in green in Table 4). The third step of the methodology consists of defining the heatwave characterisation space. The heatwave characterisation space is a space for which each dimension corresponds to one heatwave indicator. Each heatwave indicator quantifies one characteristic of the heatwaves. Thus, the third step consists of establishing a list of relevant heatwave indicators.

[7] defined a heatwave charactersation space with three indicators: the duration, the maximum temperature and the severity (also named global intensity) of the heatwave.

This three-dimensional characterisation space enables a graphical representation of a heatwave population, as shown in Figure 2. In this graph, the size of the bubbles indicates the global intensity of the heatwaves. The colour of the bubble depends on the year of the heatwave. This approach is limited: the three indicators are solely based on the evolution of the air temperature during the heatwaves. However, as mentioned in the introduction of the paper, the thermal behaviour of a building does not only result from the outside air temperature evolution. It also results from other weather parameters, such as solar radiation, sky temperature, wind speed, etc. A broader set of indicators based on those variables is needed to assess the diversity of heatwaves. Given the amount of data that is available for each heatwave, there is almost an infinite number of possible heatwave indicators. To the authors' knowledge, there was no other choice than starting the third step by the definition of an arbitrary list of indicators.

The arbitrary list of indicators proposed in the present study is composed of 3 groups.

The first group of indicators contains the heatwave duration, the day of the year corresponding to the first day of the heatwave, and some basic statistics computed from the values of the different weather variables. The basic statistics are the average value (avg), the maximum value (max), the average value of the daily maximums, the daily minima (dmin_avg/dmax_avg) and the daily amplitudes (dampl_avg), the average value during the night-time and the day-time periods (nighttime/daytime), and the average value during the week before the heatwave occurs. Here, the night-time is fixed between 10 pm and 9 am, and the daytime is fixed between 10 am and 9 pm. It is not relevant to compute all those statistics for all the weather variables. Indeed, some statistics would represent quantities that would not affect a building's thermal behaviour. Table 5 summarizes all the statistics computed for each weather variable. The The number of independent indicators measures the level of independence between the different heatwave characteristics. The higher this level of independence is, the higher the number of heatwaves that would be needed to represent the diversity of the heatwave characteristics in the dataset.

A graphical method is proposed to conduct the correlation analysis. It consists of drawing an arranged correlation matrix and identifying groups of highly correlated indicators within the matrix. This is the correlation matrix analysis (CMA). The number of indicators is then further reduced by removing insignificant indicators. These insignificant indicator values are sufficiently narrow, so that variations in their values may not lead to significant changes in building performance.

The identification of the insignificant indicators is based on the graphical representation of the distributions of indicator values. It is up to the user's degree of experience in building thermal behaviour to determine which variable may be insignificant.

Application/illustration

The values of all the indicators were computed for all the heatwaves within the extended heatwave dataset. Different combinations of linkage methods (average, ward, single, complete, weighted, centroid, median) were tested to draw the arranged correlation matrix. All those combinations resulted in similar groups of highly correlated indicators. The arranged correlation matrix on Figure 4 was obtained with the "single" linkage method corresponding to the "neared neighbour" metric between groups.

After a visual inspection of the matrix, ten groups of highly correlated indicators were arbitrarily delimited. These groups are identified with thick line rectangles in There is no real need to compute several statistics for each weather variable because those statistics are redundant. Indeed, most of the groups contain all the statistics related to one or two weather variables. For convenience, each group of indicators is named according to one of its related weather variables. The names are displayed at the bottom of Figure 4.

The heatwave duration ("duration") and the global intensity (or severity) indicators ("intensity_hwmi" and "intensity_spic") belong to the same group. This reveals that the global intensity indicators mostly measure the duration of the heatwaves and slightly modulate this duration with the air temperature. In other words, using one of the intensity indicators to describe heatwaves would not provide much additional information on the different heatwave characteristics if both the duration and the average air temperature indicators are already used. This correlation is visible in Figure 2, where it is clear that all the heatwaves with the highest intensity (larger bubbles) correspond to the longest heatwaves (on the right side of the graph).

The low correlation between solar global irradiation and solar direct irradiation was at first unexpected. It reveals the two distinct mechanisms that affect solar radiation.

Global solar irradiation is related to the period of the year at which heatwaves occur. It is at maximum when a heatwave occurs in June. It decreases when a heatwave occurs in late summer because the sun position is lower at that time. The amount of direct solar radiation seems to be more affected by the presence of clouds in the sky.

Each cooling potential indicator is based on one or two weather variables. The cooling potential indicators belong to the group of indicators related to one of their associated weather variables. The evaporative cooling potential ("cp_humid") belongs to the relative humidity group. The solar shading cooling potential ("cp_sol") is within the direct solar radiation group. The "cp_rad" indicator is highly correlated to the sky temperature group. However, it is also highly correlated with the humidity group. The correlation between the sky temperature group and the humidity group is explained by the fact that the sky vault temperature approaches the air temperature (higher values)

when the humidity content in the air is high.

For each group of indicators, one indicator was designated as the representative indicator of the group. The choice was made in an arbitrary way. The average value of the weather variable of the group was systematically chosen when this choice was made possible. The list of the independent indicators is provided in Table 6. Figure 5 shows the distributions of the ten independent indicator values. The 5 th and 95 th percentiles are given in the last two columns of Table 6. Two insignificant indicators are identified: the maximum precipitation and the main wind direction.

The 95 th percentile of the maximum precipitation values is 2.5 mm/h. According to Météo-France, precipitation can be considered light rain when its flow rate is below 3 mm/h. Thus, during heatwaves, the maximum precipitation is very low, and only very little effect of rain on building behaviour is expected.

The distribution of the main wind direction values reveals two major directions: from north to south and from south to north. These directions are opposite to each other. The wind direction mostly affects natural ventilation flow rates within cross-ventilated buildings. In most cases, the effect of the wind direction on the natural ventilation flow rate is nearly symmetrical according to the wind direction. Two opposite directions would result in a similar flow rate.

Heatwaves can therefore be characterised by a space composed of eight dimensions, which correspond to the eight significant indicators in Table 6. PCA was performed on the dataset containing the values of all the indicators for all the heatwaves. The PCA provided a list of independent linear combinations of indicators. This list can be compared to the list of independent indicators produced by the CMA.

The PCA was run until 89% of the total variance of the dataset was explained by the PCs. At this stage, eight PCs were selected. This confirms the high level of independence between the indicators. The number of PCs ( 8) is comparable to the number of independent indicators (10) selected with the CMA.

Table 7 reveals the three major contributors of the eight PCs. For almost all the PCs, the three major contributors belong to the same group of independent indicators identified with the CMA. The names of the corresponding groups are given in the last column of the table. Each group appears only once. There is a direct link between the PCs and the independent indicators. The PCA confirms the relevance of the CMA results. From the authors' point of view, the advantage of CMA, compared to PCA, is that this technique allows a visual inspection of the heatwave dataset content and selects indicators that are easier to interpret than principal components.

3.4.

Step 4: Selecting heatwaves

Objectives and method

The sampling method presented in this section aims to construct a minimal representative set of heatwaves. The minimal representative set of heatwaves covers a maximum of the diversity of the heatwave characteristics, with a minimum number of heatwaves. The sampling method consists of grouping similar heatwaves with a clustering technique and then selecting a representative heatwave for each group. Those steps are performed in the heatwave characterisation space, with normalized indicator values along the dimensions.

The grouping process is ensured by a hierarchical clustering technique. It starts with all the heatwaves being in distinct clusters. At each iteration, the two closest clusters are merged. The metric between elements for the clustering technique is the Euclidian distance in the characterization space. The linkage method is the "Ward" linkage method. This linkage method aims at producing compact, spherical clusters by selecting the clusters to be merged based on an increase in the cluster Ward variances. The two clusters that are merged are the two clusters for which the increase in the combined variance over the sum of the cluster-specific variances is the minimum.

The distance between two clusters that have just been fused is called the dissimilarity distance. At the end of the linkage process, the cumulative dissimilarity distance is displayed as a function of the number of clusters. The optimal number of clusters is a compromise between a low number of clusters and a low dissimilarity distance. It is identified graphically with the elbow method, which is commonly used with the Kmeans clustering method [START_REF] Ketchen | The application of cluster analysis in strategic management research: an analysis and critique[END_REF].

Then, a representative heatwave is selected in each cluster. To do so, the distance 𝐷 𝑘 is computed for each heatwave of the cluster. The distance 𝐷 𝑘 is the sum of the Euclidean distances between the heatwave 𝑘 and the barycentre of the other clusters:

𝐷 𝑘 = ∑ 𝑑𝑖𝑠𝑡(𝑘, 𝑏𝑎𝑟𝑦𝑐𝑒𝑛𝑡𝑒𝑟(𝑗)) 𝑗≠𝑖 ∀ 𝑘 ∈ 𝑖
The representative heatwave of a cluster is defined as the heatwave corresponding to the 85 th percentile of the distribution of 𝐷 𝑘 . This technique allows us to emphasise the dissimilarities between the different clusters and to cover a majority of the portion of the characterisation space occupied by the heatwaves without necessarily selecting outliers (Figure 6). 8 shows the number of heatwaves in each cluster. The size of the clusters is well balanced. Each cluster contains between 240 and 560 heatwaves, except the last cluster, which contains 65 heatwaves. One representative heatwave was selected within each cluster. Table 8 shows the values of the heatwave indicators for the eight representative heatwaves. A colour scale was applied to each column of the table to distinguish indicator values corresponding to lower stress (white colour) and indicator values corresponding to higher stress (red colour). The boundaries of the colour scales are the 5 th and 95 th percentiles of each indicator. Those quantiles are given in the third line of Table 8.

Each heatwave has a specific profile. The durations range between 3 days (short heatwave) and 30 days (exceptionally long heatwaves). The average air temperature oscillates between 28.8 °C and 34.3 °C. Some heatwaves have a lower average sky temperature (heatwaves 2 and 4) or a higher wind velocity value (heatwaves 2, 4 and 5), which may enhance passive cooling performance. Heatwaves 3 and 4 are the heatwaves for which the average air temperature during the preceding week is lower, enabling the building to store freshness before the beginning of the heatwave.

Heatwaves 1, 2 and 7 occur around the summer solstice, when there is plenty of solar radiation. In contrast, heatwaves 4 and 5 occur late in the summer, when there is considerably less sun. Heatwaves 3 and 6 are significantly more humid than the other heatwaves. Heatwave 5 is much drier. Finally, the sky is clear during most heatwaves, resulting in over 70% of direct solar radiation, except for heatwaves 4 and 5, where the proportion of direct radiation is 57% and 53%, respectively, and for heatwave 1, for which the direct solar radiation represents only 25% of the total solar radiation.

The eight heatwaves depict a wide variety of weather characteristics, which may result in a wide variety of building responses.

Validity of the clustering selection

The validity of the clustering selection is evaluated through a silhouette score analysis [START_REF] Rousseeuw | Silhouettes: A graphical aid to the interpretation and 625 validation of cluster analysis[END_REF], as shown in Figure 8. The silhouette score 𝑠(𝑖) of a datapoint 𝑖 lies between -1 and 1. When 𝑠(𝑖) is large (close to 1), the datapoint 𝑖 is 'well-clustered';

there is little doubt that 𝑖 has been assigned to a very appropriate cluster. When 𝑠(𝑖) is low (near zero), it is not clear whether 𝑖 should have been assigned to its cluster or to another cluster, and 𝑖 lies equally far away from both, so it can be considered an 'intermediate case'. The worst situation takes place when 𝑠(𝑖) is close to -1. In that case, it would have seemed much more natural to assign object 𝑖 to another cluster. The object 𝑖 has been 'misclassified' [START_REF] Rousseeuw | Silhouettes: A graphical aid to the interpretation and 625 validation of cluster analysis[END_REF].

The mean value of the silhouette score for the heatwave clustering is 0.06 (red dashed vertical line in Figure 8). It is quite close to zero, which indicates that clusters are overlapping.

For Clusters 1, 3, 4, 5 and 7, half of the heatwaves are 'well clustered' with a positive silhouette score, and half of the heatwaves are 'slightly misclassified', with a negative silhouette score. For Clusters 2, 6 and 8, most of the heatwaves are 'well clustered'.

Even if the silhouette score analysis reveals overlapping clusters, the representative heatwaves of each cluster seem to be quite different from each other, as represented in Table 8. It still allows a better representation of the meteorological diversity than using only one future heatwave. After performing those initial simulations, there are two ways to go deeper into the investigation of future heatwave effects on buildings. These two methods are illustrated graphically in Figure 9. The first method consists of sampling more heatwaves and performing simulations until refined statistics of the BPI values can be computed (average, standard deviation, etc.). The second usage would consist of refining the localisation of the frontier between problematic and nonproblematic heatwaves in the characterisation space. This usage would require establishing a BPI threshold to distinguish problematic and nonproblematic heatwaves.

Both ways of going further into the investigation of future heatwave effects on buildings would require a specific iterative algorithm to continue sampling heatwaves in the characterisation space.

Conclusions

The paper proposed a methodology for constructing a minimum representative set of future heatwave weather files dedicated to building performance evaluations. This methodology has been applied to the location of Lyon-St Exupery.

The set of heatwaves is representative in the sense that the heatwaves are sufficiently different from each other to cover the large spectrum of heatwave characteristics that a building could face in the location of Lyon-St Exupery. It is a minimum set; eight heatwaves were proven to be sufficient to obtain a representative set of future heatwaves.

Building simulations should be performed with the minimum representative set of heatwaves to obtain an order of magnitude of the BPI dispersion during future heatwaves and to identify to which weather characteristics the building is the most sensitive. Then, to obtain a more accurate description of the BPI value distribution or to refine the boundary between problematic and nonproblematic heatwaves, more weather files might be needed.

The minimum representative set of heatwave weather files might become an essential tool to evaluate building performance during future extreme events. It will particularly be helpful for the most robust passive cooling strategies for the future.

Figure 1 :

 1 Figure 1: Scheme of the methodology

  5 and RCP8.5, and 7 different models. The 7 models are combinations of 5 distinct GCMs and 2 distinct RCMs. Each weather projection runs from 2007 until 2097. The whole weather projection dataset represents 1260 years of weather projection.

Figure 2 :

 2 Figure 2: Heatwave characterisation space proposed by Ouzeau et al. (2016) for the heatwave detected with the "IBM" method for the RCP 4.5 scenario in the location of Lyon-Sain-Exupéry Airport. The weather predictions of all the models are mixed in this figure. The colour bar indicates the year of each heatwave. The size of the bubble represents the global intensity of the heatwave.

  second group of indicators is composed of cooling potential indicators. The cooling potential indicators quantify the opportunities offered by the weather conditions for a building to adopt passive cooling strategies. There are four "cooling potential" indicators for the four common passive bioclimatic strategies: radiative cooling (cp_rad), solar shading (cp_sol), evaporative cooling (cp_humid), and natural ventilation (cp_nv). Their definitions are represented schematically in Figure3. The radiative cooling potential is the average value of the difference between the air temperature and the sky temperature during the nighttime periods. The solar shading cooling potential is the ratio between the total horizontal direct solar radiation (rsw_dirh) and the total global horizontal solar radiation (rsw_gh). The evaporating cooling potential is the average value difference between the relative humidity of saturated air (100%) and the relative humidity of the actual outdoor air. The natural ventilation cooling potential is the average positive difference between the daily mean air temperature (DMT) and the nighttime air temperature.

Figure 3 :

 3 Figure 3: Cooling potential (by evaporation top left, by direct radiation protection bottom left, by nighttime heat radiation top right and by natural ventilation bottom right)

  An arranged correlation matrix is a correlation matrix for which the position of the indicators is arranged to put the indicators that are highly correlated closer to each other. The matrix of correlation is computed with the absolute value of the Pearson correlation. The position of the indicators in the matrix is arranged with a hierarchical clustering technique. A hierarchical clustering technique is a combination of a metric between individual elements and a linkage method. The metric between individual elements is the absolute value of Pearson's correlation. The present paper does not recommend any linkage method. The user is invited to test different methods and to select the method for which the resulting matrix allows the best visual identification of groups of highly correlated indicators. Once the groups of highly correlated indicators are identified, one indicator is retained within each group. The other indicators are removed from the list. They are redundant.

Figure 4 .

 4 Figure 4. Two indicators are not correlated with each other and thus form two groups of one indicator: wind direction and average air temperature one week before the heatwave.

Figure 4 :

 4 Figure 4: Correlation matrix (in absolute values) between the different heatwave indicators. Strong correlations are shown in red, and weak ones in blue. Black rectangles represent the groups of highly correlated indicators. The names of the groups are given at the bottom of the figure. Dashed rectangles represent insignificant groups.

Figure 5 :

 5 Figure 5: Histograms of independent indicator values among the extended heatwave dataset

Figure 6 :

 6 Figure 6: Representative heatwave selection and the space it covers on a 2D projection of the characterisation space.

3. 4 . 2 .

 42 Application/illustrationThe result of the clustering process for Lyon-St Exupery is represented as a dendrogram in Figure7(left graph). In this graph, each junction between two branches represents the fusion between two clusters. The vertical distance between two consecutive junctions is the dissimilarity distance.

Figure 7 (

 7 Figure7(right graph) shows the evolution of the dissimilarity distance as a function of the number of clusters. The curve is strictly decreasing: the higher the number of clusters is, the more similar the heatwaves are in each cluster. The elbow of the curve occurs when the number of clusters is equal to eight. Onwards, the decrease in the dissimilarity distance is strongly reduced. Choosing a larger number of clusters would not significantly reduce the distance between clusters, but it would considerably increase the number of simulations for assessing building performances. The last column of Table8shows the number of heatwaves in each cluster. The size of the

Figure 7 :

 7 Figure 7: Dendrogram of the heatwave dataset in the characterisation space (left graph). Euclidean distance as a function of the number of clusters (right graph). The red dashed lines show the selected number of clusters and the associated Euclidean distance.

Figure 8 :

 8 Figure 8: Silhouette score for each heatwave in the dataset. The red vertical dashed line indicates the mean of the silhouette score. The number on the right of each cluster indicates the mean value of the silhouette score for each cluster.

4 .

 4 Usage and beyond530The minimum representative set of heatwaves should be the starting point of the evaluation of building performances under future heatwaves. Since the minimum representative set of heatwaves covers the whole diversity of future potential heatwave characteristics, the simulations performed with this set of weather files may reproduce the expected variety of building thermal behaviour during future heatwaves. Those simulations may also produce a first population of building performance indicator values (BPI) and allow us to have a first estimation of the amplitude of those values.

Figure 9 :

 9 Figure 9: Dispersion of a building performance indicator due to the variability of heatwaves (left graph) and frontier between problematic and nonproblematic heatwaves in a 2-dimensional projection of the characterisation space (right graph).

Table 1 :

 1 General information on the case study

	Name	Latitude	Longitude	Altitude	Time Zone
	Lyon-Saint Exupéry	45.7261	5.09083	250 m	UTC +1:00 (Europe/Paris)
		45° 43' 34'' North	5° 5' 27'' East		Summer time: UTC +2:00
					Winter time: UTC +1:00
	3. Detailed description of the methodology
	3.1. Step 1: Collecting climate projection data
	3.1.1. Objectives and method

Table 2 :

 2 Minimum requirements for a CORDEX weather projection to be considered suitable

	List of required variables	Maximum	Bias-correction availability
		Timestep
	Dry air temperature [K]	3 hrs	Yes
	Specific humidity [kg/kg]	3 hrs	-
	Atmospheric pressure [Pa]	3 hrs	-
	Cloud fraction [%]	3 hrs	-
	Global horizontal radiation [W/m²]	3 hrs	Yes
	Wind speed [m/s]	3 hrs	Yes

Table 3 :

 3

	GCM	RCM	Ensemble	rcp45	rcp85

List of the weather projections collected from the CORDEX database. Crosses indicate data availability.

Table 4 :

 4 Average annual number, air temperature and duration of heatwaves detected from each detection method.

	Detection	Average annual number of	Average air	Average	Total number of
	methods	heatwaves	temperature	duration	heatwaves
	Spic, Sdeb, Sint	1.89	28.3 °C	14.5 d	2384
	IBMn, IBMx	1.77	30.3 °C	6.3 d	2229
	EHF	6.06	24.8 °C	11.1 d	7641
	HWMI	6.65	25.6 °C	6.9 d	8375
	3.3. Step 3: Heatwave characterisation

3.3.1. Objectives and method

Table 5 :

 5 First group of the extended set of indicators, inventory of the statistics computed for each weather variable

	Weather variable	avg	max	dampl_avg	dmin_avg/dmax_avg	nighttime/daytime
	Dry air temperature [°C]	x	x	x	x	x
	Specific humidity [kg/kg]	x
	Relative Humidity [%]	x				x
	Sky temperature [°C]	x				x
	Global horizontal radiation [W/m²]	x
	Direct horizontal radiation [W/m²]	x
	Cloud cover	x
	Wind speed [m/s]	x				x
	Precipitation [kg/m²/s]	x	x

Table 6 :

 6 Minimum and maximum values of each independent indicator.

	Indicator	Variable name	Unity	Q5	Q95	Significant
	Maximum precipitation	pr_max	mm/h	0	2.5
	Average air temperature	ta_avg	°C	28.36	32.99	X

Table 7 :

 7 Major contributors of the PCs

		1 st	2 nd	3 rd	Clustering group
	PC1	ta_daytime	ta_avg	ta_dmax_avg	Air temperature
	PC2	tsky_avg	tsky_nighttime	tsky_daytime	Sky temperature
	PC3	ws_avg	ws_nighttime	ws_daytime	Wind
	PC4	rsw_dirh_avg	day_of_year	rsw_gh_avg	Solar radiation
	PC5	ta_dampl_avg	rsw_gh_avg	day_of_year	-
	PC6	duration	intensity_hwmi	intensity_spic	Duration and intensities
	PC7	pr_avg	pr_max	cp_sol	Rain
	PC8	ta_avg_prior	wd_main	day_of_year	-

Table 8 :

 8 Characteristics of the heatwaves contained in the minimum representative set of heatwaves

		duration ta_avg	tsky_avg	ws_avg ta_avg_prior rsw_gh_avg rh_avg	cp_sol
	Unit	[𝒅]	[°𝑪]	[°𝑪]	[𝒎/𝒔]	[°𝑪]	[𝑾/𝒎²]	[%]	[-]	nb heatwaves
	range [3; 16] [28.4; 33] [10.1; 18.6] [1.82; 4.37] [22.2; 30.4] [207; 321] [27; 57.3] [0.4; 0.9]
	1	3	31.16	17.68	3.65	25.18	307.88	41.4	0.25	250
	2	4	30.89	13.7	5.46	24.62	299.19	31.78	0.72	246
	3	5	29.00	17.35	1.90	23.14	280.1	59.76	0.88	297
	4	6	28.79	10.1	3.96	23.14	236.33	41.16	0.57	564
	5	6	33.73	14.62	3.96	26.45	220.03	22.6	0.53	253
	6	6	29.31	15.43	2.24	23.75	284.94	53.3	0.91	238
	7	10	34.31	17.17	3.42	24.82	300.13	27.24	0.8	316
	8	30	32.81	18.57	2.8	31.57	265.86	36.71	0.8	65

Acknowledgements

We acknowledge the World Climate Research Program's Working Group on Regional Climate and the Working Group on Coupled Modelling, the former coordinating body of CORDEX and the responsible panel for CMIP5. We also thank the climate modelling groups (listed in Table 3 of this paper) for producing and making available their model output. We also acknowledge the Earth System Grid Federation infrastructure, an international effort led by the U.S. Department of Energy's Program for Climate Model Diagnosis and Intercomparison, the European Network for Earth System Modelling and other partners in the Global Organisation for Earth System Science Portals (GO-ESSP). The authors gratefully acknowledge the Lyon Urban School (LUS) for funding this research work (ANR-17-CONV-0004).


Yosra Gargouri
email: yosra.gargouri@cea.fr

Nicolas Ravot
email: nicolas.ravot@cea.fr

Mariem Slimani
email: mariem.slimani@cea.fr

Mickael Cartron
email: mickael.cartron@cea.fr

Enhancing Reflectometry Systems with CHIRP-OMTDR and Compressed Sensing: A Study on Signal Recovery Quality

Keywords: OMTDR, CHIRP-OMTDR, Random Demodulator, Reflectometry, Compressed Sensing

This paper investigates the application of Compressed Sensing (CS) in Orthogonal Multi-tone Time-Domain Reflectometry (OMTDR) using a specific signal known as CHIRP-OMTDR. The Random Demodulator (RD) has been chosen as the Analog-to-Information Converter (AIC) in the acquisition chain. The study focuses on the influence of RD's filter selection, signal length, and compression factors on signal recovery quality. The results demonstrate that reconstruction quality improves with longer signal lengths, and that the choice of low-pass filter type and order has minimal impact. The findings contribute to the enhancement of reflectometry systems based on Compressed Sensing.

I. INTRODUCTION

Reflectometry is a nondestructive sensing system used to identify and characterize electrical faults in cables or transmission lines. It involves transmitting a signal through a medium and analyzing the reflections caused by impedance changes. Reflectometry has also other applications such as detecting damage in photovoltaic systems, characterizing nanostructures, and evaluating moisture content in soils [START_REF] Kingston | A sstdr methodology, implementations, and challenges[END_REF].

A significant challenge in reflectometry is enhancing the accuracy of locating and characterizing impedance discontinuities, which contain valuable information. This task requires the use of high-frequency Analog-to-Digital Converter (ADC). However, the cost, power consumption, and complexity associated with very high-frequency ADCs introduce significant limitations. As a result, alternative solutions like compressed sensing have emerged to overcome these challenges [START_REF] Tzila Ajamian | A novel compressive sampling approach for detecting hard defects in complex wire networks[END_REF].

Compressed sensing theory enables accurate signal recovery with fewer samples than what the Nyquist-Shannon theorem dictates. This reduction in required samples allows for a decrease in sampling frequency and memory requirements. Previous studies [START_REF] Tzila Ajamian | A novel compressive sampling approach for detecting hard defects in complex wire networks[END_REF] have demonstrated promising results in detecting impedance discontinuities using compressed sensing on linear chirp signals. This study goes a step further by investigating OMTDR signals and examining the impact of filter selection, signal length, and compression factors on the quality of the reflected signal reconstruction. The findings will contribute to a better understanding of the factors influencing signal reconstruction, ultimately improving the accuracy and reliability of CS-based reflectometry systems.

II. THEORY OF COMPRESSED SENSING A. Principle of Compressed Sensing

In traditional Digital Signal Processing, signals are typically sampled uniformly following the Nyquist-Shannon theorem, where the sampling rate is at least twice the signal bandwidth. Compressed Sensing, also known as Compressive Sensing, Compressive Sampling or CS, is a data acquisition technique that aims to capture signals directly in a compressed form with nonadaptive measurements.

The main idea behind CS is that many real-world signals exhibit sparsity, meaning that they can be represented by a sparse or compressible set of coefficients in a suitable basis. By leveraging the knowledge of signal sparsity, the compressed measurements obtained through CS can be efficiently reconstructed using advanced signal processing algorithms. It becomes possible to accurately reconstruct the original signal at a significantly lower rate than the Nyquist rate. This reduction in sampling rate has various advantages, including reduced storage requirements, data transmission bandwidth, and power consumption during acquisition [START_REF] Gargouri | Compressive sampling for efficient astrophysical signals digitizing: From compressibility study to data recovery[END_REF].

B. Analog-to-Information Converter (AIC)

The concept of compressed sensing has led to the development of a new type of converter known as the Analog-to-Information Converter (AIC) [START_REF] Rapuano | Analog-to-information converters: research trends and open problems[END_REF]. AICs differ from conventional ADCs in that they can operate at lower sampling rates than those dictated by the Nyquist-Shannon theorem thanks to the CS techniques. Prominent AIC architectures include the Non-Uniform Sampler (NUS), Random Demodulator (RD), and Random Modulator Pre-Integrator (RMPI) [START_REF] Gargouri | Compressive sampling for efficient astrophysical signals digitizing: From compressibility study to data recovery[END_REF].

III. CS-BASED REFLECTOMETRY ARCHITECTURE A. Conventional reflectometry vs CS-based reflectometry

The conventional reflectometry chain is depicted in Fig. 1. A signal is injected into the device under test (DUT) using a digital-to-analog converter (DAC). This injected signal propagates through the medium and reflects upon encountering impedance discontinuities. The reflected signal is then acquired using an analog-to-digital converter (ADC) operating at a sampling frequency f s . Following the acquisition, signal processing techniques are employed for fault detection and localization purposes.

By utilizing CS, the ADC operating at f s in the conventional chain will be replaced by an Analog-to-Information Converter (AIC) operating at a lower frequency, f AIC < f s . This change necessitates a signal reconstruction phase, which is followed by signal processing for fault detection and location, as shown in Fig. 2.

In this study, we will focus on exploring the feasibility of implementing Compressed Sensing technology specifically for OMTDR method.

B. OMTDR and CHIRP-OMTDR

Reflectometry-based techniques are categorized based on the injected waveforms and the methods used for analyzing the reflected signal. The analysis of the reflected signal can be performed in either the time domain (Time Domain Reflectometry) or the frequency domain (Frequency Domain Reflectometry). Various well-known reflectometry methods are reviewed in [START_REF] Cynthia M Furse | Fault diagnosis for electrical systems and power networks: A review[END_REF].

Previous studies on CS-based reflectometry mainly focused on Time-Domain Reflectometry (TDR) methods using linear chirp signals [START_REF] Tzila Ajamian | A novel compressive sampling approach for detecting hard defects in complex wire networks[END_REF]. In this paper, we will explore the Orthogonal Multi-Tone Time Domain Reflectometry (OMTDR) method.

OMTDR offers not only real-time diagnosis capabilities but also facilitates data transmission in complex wiring networks [START_REF] Ben Hassen | Smart connectors fusion for incipient defect diagnosis on aging aircraft wiring systems[END_REF]. In OMTDR, a signal consisting of multiple tones or frequencies is transmitted, with these tones being encoded using the M-phase shift keying (M-PSK) digital modulation technique to carry information [START_REF] Cynthia M Furse | Fault diagnosis for electrical systems and power networks: A review[END_REF]. However, the classical OMTDR signal is not sparse in the time domain (Fig. 3(a)), frequency domain, or time-frequency domain (Fig. 3(b)). This poses a challenge since CS requires signals to exhibit sparsity in a known domain.

To address this challenge, we have developed a specific OMTDR signal, referred to as CHIRP-OMTDR (Fig. 3(c)), which follows the M-PSK modulation and features a compressible chirp waveform in the time-frequency domain (Fig. 3(d)). Detailed information on generating this signal can be found in patent application [START_REF] Yosra | Method for generating a multi-carrier reflectometry signal having the form of a" chirp" signal[END_REF]. The CHIRP-OMTDR signal format enables us to explore the application of Compressed Sening in OMTDR-based reflectometry systems and delve deeper into its potential benefits.

IV. AIC ARCHITECTURE MODELING

A. Choice of the AIC architecture

We have chosen the Random Demodulator as the AIC architecture in the acquisition chain. The choice of this architecture is motivated by its efficiency on compressible signals in the time-frequency domain for linear chirp waves [START_REF] Tzila Ajamian | A novel compressive sampling approach for detecting hard defects in complex wire networks[END_REF] and its relative simplicity of implementation compared to an RMPI architecture.

The RD architecture, depicted in Fig. 4, comprises three main components: a mixer, a filter, and an ADC.

In the RD-based AIC, the reflected signal (S ref ) is demodulated at the device's input using a pseudo-random sequence (PN) composed of +1 and -1 values. This demodulation process operates at a frequency that is equal to or higher than the Nyquist frequency. Subsequently, the demodulated signal is filtered and then sub-sampled using an ADC clocked at a frequency f AIC , where f AIC equals f s /CF (CF denotes the compression factor employed).

B. Choice of filter and compression factor

In the literature, various types of filters have been employed in Random Demodulator. Historically, early filters utilized integrator designs [START_REF] Laska | Theory and implementation of an analog-to-information converter using random demodulation[END_REF] [START_REF] Ragheb | A prototype hardware for random demodulation based compressive analog-to-digital conversion[END_REF]. Butterworth filters have also been utilized, including second-order Butterworth low-pass filters [START_REF] Jia | Design and analysis of a hardware-efficient compressed sensing architecture for data compression in power quality data acquisition[END_REF] and fourth-order Butterworth low-pass filters [START_REF] Pawel | Sensitivity of the random demodulation framework to filter tolerances[END_REF].

In order to determine the appropriate filter for our architecture, we simulate different filters including the integrator and Butterworth filters of orders one to six. We compare their reconstruction accuracy based on the Signal-to-Error Ratio (SER). Additionally, we investigate the impact of signal lengths, and compression factor on the system's performance. To achieve this, we proceed as follows:

1) During the emission phase, we generate the CHIRP-OMTDR signal with a variable length (N, ranging from 128 to 2048 samples) at a frequency of 200 MHz.

2) The signal is injected and propagated into a coaxial cable.

3) The reflected signal passes through the AIC block: it is multiplied by a pseudo-random sequence at 200 MHz. 4) Next, it undergoes filtering. We explore different filter choices, including an integrator filter and Butterworth filters ranging from the first to the sixth order. 5) Subsequently, the filtered signal is downsampled using compression factors of 2, 4, and 8. 6) Finally, we reconstruct the signal and evaluate the reconstruction accuracy (SER) using the general formula. SER = 20 • log 10

∥R CV ∥2 ∥(R CV -R CS )∥2
where R CV represents the reflectogram of a conventional reflectometry chain and R CS represents the reflectogram of a CS-based reflectometry chain. Figures 5,6 and 7 illustrate the variation of SER as a function of the signal length N and the nature of the filter for compression factors 2, 4 and 8, respectively.

We observe that irrespective of the filter type and its order, there is a consistent trend in terms of reconstruction quality, with the results being very similar. Therefore, what matters most from our standpoint is not necessarily having a high-fidelity low-pass filter, but rather having a filter with a stable impulse response. This is crucial due to the sensitivity of the Random Demodulation (RD) technique to distortion introduced by the low-pass filter component [START_REF] Pawel | Sensitivity of the random demodulation framework to filter tolerances[END_REF].

Furthermore, we observe that, for the same compression factor CF, the reconstruction quality (SER) improves as the signal length N increases. This finding aligns with the state of the art in CS, which suggests that the minimum measurement requirement for accurate reconstruction scales approximately logarithmically with the signal size. In other words, as the signal length increases, the number of required measurements for accurate reconstruction also increases, but at a slower rate than a linear relationship. However, it is worth noting that the reconstruction quality deteriorates when using a compression factor (CF) of 8. We will exclude this compression factor and instead focus on CF = 2 and CF = 4.

C. Choice of the length of the injection signal

To determine the optimal signal length N, two factors must be considered: signal reconstruction error precision (SER) and the time required for signal reconstruction. To determine the optimal signal length, we conducted a statistical analysis of 30 simulations using two filter types (an integrator and a 6thorder Butterworth) and two compression factors (CF=2 and CF=4). By applying Compressed Sensing with varying signal lengths (ranging from 128 to 2048), we reconstructed the signal, measured the reconstruction time, and evaluated the SER.

Our observations revealed that increasing the signal length improves signal reproduction accuracy, resulting in a higher SER. However, longer sequences also lead to longer execution times during the reconstruction phase. To strike a balance between SER and reconstruction time, we found that selecting N=512 points provides a well-balanced compromise.

Fig. 9 illustrates an example of the reconstructed reflectogram quality with a reconstruction precision of 8.42 dB compared to a standard reflectogram without CS. This reconstruction precision was evaluated after applying CS to a 512-point signal, using a 6th-order Butterworth filter and subsampling by a compression factor of 4, as compared to the standard case.

V. CONCLUSION

In conclusion, our study has demonstrated the effectiveness of a CS-based reflectometry architecture in accurately recovering reflected signals in the context of OMTDR. By employing a compressed factor of 4 and a 512-point CHIRP-OMTDR signal, we achieved reliable signal recovery. Furthermore, our findings indicate that the selection of low-pass filter type and order has negligible influence on the reconstruction quality. After completing the phase of simulations and testing, our next step involves conducting an experimental test bench using evaluation boards. The ultimate goal is to develop an analog reflectometry printed circuit board (PCB) based on Compressed Sensing.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Schematic diagram of conventional reflectometry architecture

Fig. 3 .Fig. 4 .

 34 Fig. 3. (a) The time representation of the classical OMTDR signal. (b) the time-frequency representation of the classical OMTDR signal. (c) the time representation of the OMTDR signal in chirp form (d) the time frequency representation of the OMTDR signal in chirp form.

Fig. 5 .

 5 Fig. 5. Evaluation of the SER as a function of the number of signal points and the nature of the filter (cutoff frequency fc=50 MHz for Butterworth filters) for a compression factor of 2.

Fig. 6 .

 6 Fig. 6. Evaluation of the SER as a function of the number of signal points and the nature of the filter (fc=25 MHz for Butterworth filters) for a compression factor of 4.

Fig. 7 .

 7 Fig. 7. Evaluation of the SER as a function of the number of signal points and the nature of the filter (fc=12.5 MHz for Butterworth filters) for a compression factor of 8.

Fig. 8 .

 8 Fig. 8. Signal Length Estimation.

Fig. 9 .

 9 Fig. 9. Reflectogram generated with the conventional reflectometry architecture at 200 MHz (blue curve), reflectogram generated with the reflectometry architecture based on CS at 50 MHz (N=512 samples, CF=4, Butterworth filter of order 6 and fc=25 MHz).

ACKNOWLEDGMENT

The authors gratefully acknowledge funding support from BPI France and his partners Nicomatic and Heatself for the project SECURE-HOP.


Ngac Wenjing Zhang

Ky Nguyen

Eric Semail

Yanliang Xu

Axial Stress Analysis and Comparison of the Novel Dual 3-phase Axial Flux Permanent Magnet Machines

come

Abstract--In previous research, a novel three-phase dualstator axial flux permanent magnet machine characterized by the advantages of compact structure and low moment of inertia is proposed for industrial robot application. In order to improve its functional reliability furtherly, the dual threephase axial flux permanent magnet machine (DTP-AFPM) is firstly proposed. Benefiting from the combination of 12 slots/ 10 poles, the coil configuration of one stator disk can be modified to a dual three phase full-pitch winding straightforwardly, as a result, one kind of the DTP-AFPMs is achieved which is named as the no shift model in this paper. For eliminating the coil reconfiguration on each stator disk and the connection of the coils belonging to the same phase between two stator disks, the shift model which is based on the shift of the two stator disks to obtain the phasor difference between two three-phase windings is introduced. The characteristics and electromagnetic performances of these two models are analyzed and compared. However, it should be noted that the light-weight disk-type rotor of DTP-AFPMs also degrade the strength of the rotor. The proposed DTP-AFPMs are more sensitive to the axial stress on the rotor which introduces not only the vibration and noise but also the deformation or even the damage of the rotor. Thus, the axial stress on the rotor is investigated and treated as a critical evaluation indicator. The axial stress is analyzed under both healthy and fault conditions and its distribution on the rotor is given on a 2-D plane.

Index Terms--Axial flux permanent magnet machine, dual three-phase, one-phase-open, fault tolerant control, axial stress.

I. INTRODUCTION

In previous research [START_REF] Zhang | Research on a Novel Transverse Flux Permanent Magnet Motor With Hybrid Stator Core and Disk-Type Rotor for Industrial Robot Applications[END_REF], a novel axial flux permanent magnet machine (AFPM) as shown in Fig. 1 is proposed for the industrial robot application where the low volume, low moment of inertia and additional space for kinds of accessories are required.

Normally, the 3-phase machine is enough for most industry applications. However, due to its unsatisfactory presentation on reliability, the demands for the multiphase machine has been growing. Besides the superiority in reliability [START_REF] Liu | A novel multiphase brushless powersplit transmission system for wind power generation[END_REF][3], the multiphase machine also combines the advantages of low voltage rating [START_REF] Levi | Multiphase electric machines for variable-speed applications[END_REF], high power density [START_REF] Zhang | Design and comparison of three-phase and five-phase FTFSCW-IPM motor openend winding drive systems for electric vehicles applications[END_REF] and high torque density by injecting harmonics [START_REF] Wang | Dual-rotor multiphase permanent magnet machine with harmonic injection to enhance torque density[END_REF]. In [START_REF] Locment | Vectorial approachbased control of a seven-phase axial flux machine designed for fault operation[END_REF], a seven-phase axial flux permanent magnet machine supplied by a seven-leg voltage source inverter is proposed, the machine is able to work easily with one or two open-circuited phases by means of two supplementary degrees of freedom for fault mode configuration. In [START_REF] Slunjski | General torque enhancement approach for a nine-phase surface PMSM with built-in fault tolerance[END_REF], the authors analyze a nine-phase surface-mounted PMSM, the output torque is improved by using non-sinusoidal back EMF and the sensor fault can be solved based on the harmonic components in EMF. In [START_REF] Zhao | Multiphase PMSM with asymmetric windings for more electric aircraft[END_REF] and [START_REF] Villani | High reliability permanent magnet brushless motor drive for aircraft application[END_REF], the authors both develop multiphase permanent magnet machine for aircraft application considering the strict requirement for the safety of the flight which can be improved by the faulttolerant capability of multiphase winding structure. Dual 3-phase machine is one of the widely applied form which combines the advantages of control easily for conventional 3-phase machine and fault-tolerant capability for multiphase machine. Besides, the proposed novel AFPM is more suitable to be transformed to a dual 3-phase machine due to its winding configuration and dual-sides stator structure [START_REF] Zhang | Research on a Novel Transverse Flux Permanent Magnet Motor With Hybrid Stator Core and Disk-Type Rotor for Industrial Robot Applications[END_REF].

The AFPM benefits from the advantages of compact structure [START_REF] Shokri | Comparison of performance characteristics of axial-flux permanent-magnet synchronous machine with different magnet shapes[END_REF] and high torque density [START_REF] Simón-Sempere | Optimisation of Magnet Shape for Cogging Torque Reduction in Axial-Flux Permanent-Magnet Motors[END_REF]. However, the AFPMs is generally characterized by a disk-type rotor of which the strength is relatively lower compared with the radial flux rotor. Thus, the analysis of the axial stress on the rotor is important. According to the Maxwell stress tensor theory, there are three types of stress exerted on the rotors of the AFPM: tangential stress which is the origin of the output torque, radial stress and axial stress that are the results of the vibration and acoustic noise. And it is proved that the magnitude of the radial stress in AFPM is much lower compared with the axial stress and can be neglected as illustrated in [START_REF] Park | A numerical prediction model for vibration and noise of axial flux motors[END_REF]. In terms of the single-side AFPM, the axial force on rotor is dangerously high. And the uneven axial force distribution on the surface of the rotor also deteriorate the problems of vibration and noise [START_REF] Deng | Axial force and vibroacoustic analysis of external-rotor axial-flux motors[END_REF]. To deal with this problem, the double-stator single-rotor symmetric topology is adopted of which the axial force on the rotor can be balanced. However, the balance will be destroyed while there is an asymmetrical condition either in mechanical structure or the winding distribution and it is bound to deteriorate furtherly while the AFPM are working under the post-fault condition [START_REF] Lu | Reduction of unbalanced axial magnetic force in postfault operation of a novel six-phase double-stator axial-flux PM machine using model predictive control[END_REF]. For the proposed structure in [START_REF] Zhang | Research on a Novel Transverse Flux Permanent Magnet Motor With Hybrid Stator Core and Disk-Type Rotor for Industrial Robot Applications[END_REF] which adopt a coreless lightweight rotor, this problem is more significant. As a result, the axial stress on the rotor is a vital indicator to evaluate the feasibility of the proposed novel structure.

In this paper, the transformation from a previously proposed 3-phase AFPM to two kinds of dual 3-phase axial flux permanent magnet machines (DTP-AFPM) is introduced firstly. The electromagnetic performances are analyzed. In section III, the control method under both healthy and fault conditions are presented briefly. One phase open-circuit fault is considered and the fault tolerant control (FTC) strategy is adopted. The simulation results are given and the phase current is used to calculate the axial stress. In section IV, the axial stress distribution on the rotor is analyzed and compared between the two DTP-AFPM models in three different stages.

II. STRUCTURE AND ELECTROMAGNETIC ANALYSIS

A. Structure and Coil Configuration of DTP-AFPMs

The structure of 3-phase AFPM that has been researched in previous work is shown in Fig. 2. The double-stator single rotor topology is adopted of which the two stator disks are symmetrically arranged at the two sides of the rotor disk as shown in Fig. 2(a). Each stator disk is composed of 12 segmented armatures. As shown in Fig. 2(b), the segmented armature consists of a stator shoe made by SMC material, stator body made by Si-steel and the coil. The commonly used coil configuration W1 for 10 poles/12 slots machine as shown in Fig. 3(a) is used by which a short-pitch concentrated winding is obtained. The coil configurations of both two stator disks are identical. The phasor arrangement of the two stator disks is shown in Fig. 4(a). All coils of two stator disks belong to the identical phase are connected serially.

Based on the 10 poles/12 slots combination of 3-phase AFPM, it is convenient to turn it into a novel DTP-AFPM which benefits from some additional advantages such as fault-tolerance capability, lower voltage-level and lower torque ripple. A dual 3-phase full-pitch concentrated winding is achieved by modifying the coil configuration in accordance with W2 as given in Fig. 3(b). It can be noticed that the phasor difference of two winding sets is 30 electrical degrees which eliminates the 6th harmonics in torque ripple. The phasor arrangement of the two stator disks is shown in Fig. 4(b). The coils of the identical phase sequence are connected serially. Based on the 10 poles/12 slots combination of 3-phase AFPM, it is convenient to turn it into a novel DTP-AFPM which benefits from some additional advantages such as fault-tolerance capability, lower voltage-level and lower torque ripple. A dual 3-phase full-pitch concentrated winding is achieved by modifying the coil configuration in accordance with W2 as given in Fig. 3(b). It can be noticed that the phasor difference of two winding sets is 30 electrical degrees which eliminates the 6th harmonics in torque ripple. The phasor arrangement of the two stator disks is shown in Fig. 4(b). The coils of the identical phase sequence are connected serially.

Besides, a more convenient method to obtain a dual 3phase AFPM with the elimination of the coil reconfiguration and the electrical connection between two stator disks is proposed. The coil configurations of both stator disks continue to adopt W1 while the back stator is shifted by certain mechanical degrees as shown in Fig. 5. Then, two sets of 3-phase windings with phasor difference are obtained on the two stator disks as shown in Fig. 4(c). The shifted angle is set to 6 mechanical degrees (30 electrical degrees).

B. Comparison of the Electromagnetic Analysis Results

According to the design parameters given in [START_REF] Zhang | Research on a Novel Transverse Flux Permanent Magnet Motor With Hybrid Stator Core and Disk-Type Rotor for Industrial Robot Applications[END_REF] and the modification illustrated above, two DTP-AFPM models (shift model and no shift model) are constructed and analyzed using ANSYS MAXWELL.

The back-EMF waveforms of these two models are given in Fig. 6. The models are operated as generators and the open-circuit output phase voltages are calculated. For the shift model, the RMS of phase voltage is 79.5V. For no shift model, it is 82.4V. The harmonic component in back-EMF of the shift model and no shift model is analyzed using Fourier analysis and compared as shown in Fig. 7. The total harmonic distortion (THD) of the phase back-EMF of the shift model is 8.4% while it is 11.1% for the no shift model. It can be concluded that the no shift model is with a higher back-EMF benefiting from which it can produce higher torque under the same input current, however, its harmonic component is more significant. The main reason of the difference on the back-EMF between these two models is that the winding factor of the no shift model is 1 that means there is no discount for the synthesis of back-EMF and no offset of the specific harmonic components. The main harmonic component in back-EMF is 3 rd as shown in Fig. 7 which can be eliminated eventually due to the separated neutral points. The rated torque of DTP-AFPMs is 17.2Nm. The comparison of the percentage of the cogging torque relative to rated torque between the shift and no shift models is shown in Fig. 8. The percentage value of the shift model is while it is of the no shift model.

III. CONTROL OF DTP-AFPM

The control system of DTP-AFPM is shown in Fig. 9. The mathematical model of DTP-AFPM is constructed based on the vector space decomposition (VSD) method in [START_REF] Kallio | Decoupled dq model of double-star interior-permanent-magnet synchronous machines[END_REF]. By use of the VSD, all the variables of the machine are transformed into a decoupled frame consists of the orthogonal D1-Q1 axes and D2-Q2 axes. The torque equation of DTP-AFPM under decoupled DQ-frame is simplified in the form of ( 2) and the current components on D1, D2 as well as Q2 axis are equal to 0 while the machine operates under healthy condition and the harmonics in phase current is neglected. While the DTP-AFPM operates under a fault condition in which the open-circuit fault in phase A2 is assumed, the output electromagnetic torque will fluctuate greatly if the reference current is still calculated using the control strategy under healthy condition. In order to output a smooth torque under fault condition, a fault tolerant control (FTC) strategy in [START_REF] Feng | Open-phase fault modeling and optimized fault-tolerant control of dual three-phase permanent magnet synchronous machines[END_REF] is referenced. The reference current on DQ-frame under fault condition is reconfigured as given in (3). The factor F I is determined by the torque equation ( 4).

F D f Q f F D f F F Q f I i i I i I I i π θ π θ π θ  = -   =    = +    = -   , (3)
in which symbol θ indicates the rotor electrical position.

1, [START_REF] Liu | Design and analysis of a new fault-tolerant permanent-magnet vernier machine for electric vehicles[END_REF] 3

em M Q f M F T p i p I Ψ Ψ = = (

IV. ANALYSIS OF AXIAL STRESS

As an axial flux machine characterized by the lightweight disk-type coreless rotor, DTP-AFPM is more sensitive to the axial stress. To prevent the potential possibility of unacceptable NVH (Noise, Vibration and Harshness) and mechanical damage of the machine during the operation, the axial stress on the rotor is significant to be analyzed and treated as a critical evaluation indicator for judging the feasibility of the proposed DTP-AFPM. It is clear that the axial stress on the rotor of the no shift model is balanced while it is unbalanced in the shift model due to its asymmetrical structure. In this section, the axial stress and total axial force on the rotor are calculated and compared for both shift model and no shift model.

A. Calculation Method of Synthetic axial stress

According to the Maxwell Stress Tensor (MST) principle, the normal stress on one face can be calculated by following equation: The detailed calculating procedure of the synthetic axial stress under both healthy and fault condition is as follows:

2 2 2 0 1 ( ) 2 z z x y f B B B µ = -- , (5)
1) Calculate the phase current of DTP-AFPM under different conditions; 2) Impose the calculated phase current as shown in Among the six subgraphs in Fig. 15, there is no obvious regularity can be found. The peak value of each figure is indicated and the extremum is about 0.1MPa which is most likely due to the calculation errors.

C. Synthetic Axial Stress in the Shift Model

The rotor axial forces of the shift model while the machine operates in different stages are shown in Fig. 16. The average value of the rotor axial force in one electrical circle is -2.24N under healthy condition (stage 1). It means that the bearing on one side will support the axial force throughout. Under fault condition (stage 2), the rotor axial force deteriorates sharply that fluctuates between ± 100N. It is bound to cause the serious vibration and noise. Under fault condition with FTC (stage 3), the fluctuation of the rotor axial force waveform decreases clearly. However, it is still at a higher level.

In Fig. 16, a special rotor displacement position =294° is highlighted where the axial force on the rotor reaches the extremum under fault conditions (stage 2 and 3). The synthetic axial stress distribution of the shift model at =294° is shown in Fig.

V. CONCLUSION

The DTP-AFPM is proposed in this paper including the shift model and the no shift model. The electromagnetic performances of these two models are compared. The synthetic axial stress on the rotor of the DTP-AFPM is analyzed. The main conclusions according to the investigations in this paper are summarized as follows:

1) The shift model is a more accessible structure for DTP-AFPM due to the elimination of both the electrical connection between the two stator disks and reconfiguration of the coils.

2) The no shift model is characterized by a higher back-EMF benefiting from its full-pitch winding arrangement, while its harmonics in back-EMF is degraded compared with the shift model. The cogging torques of the two models are with the same magnitude and reversed phasor.

3) The fluctuation of the total axial force on the rotor of the shift model is much higher than that of the no shift model, especially operating under the fault condition. Although the FTC strategy can greatly reduce the axial force on the rotor of the shift model, its influence on the vibration and noise still exists.

4) The breaking possibility of the PMs of the shift model is much higher according to the synthetic axial stress distribution diagram given in section VI.

5) Considering the axial stress is the crucial indicator for the mechanical safety of the proposed DTP-AFPM, the no shift model is the best choice.

Fig. 1 .

 1 Topology of the proposed novel AFPM in [1]: (a) Section view (b) Explosive view.

2 .

 2 Structure of previous 3-phase DTFPM: (a) Total structure (b) Single armature stator core (c) Rotor part. (a) (b) Fig. 3. Coil configuration of one stator disk: (a) W1: 3-phase shortpitch concentrated winding (b) W2: dual 3-phase full-pitch concentrated winding. (The digital indicates the number of the stator core as shown in Fig. 2)

Fig. 4 .Fig. 5 .

 45 Fig. 5. Physical diagram of the stator disks with shift.

Fig. 6 .

 6 Fig. 6. Comparison of back-EMF waveforms.

Fig. 7 .

 7 Fig. 7. Comparison of the harmonic component in back-EMF.

Fig. 8 .

 8 Fig. 8. The comparison of the percentage of the cogging torque relative to rated torque.

Fig. 9 .

 9 Fig. 9. The schematic diagram of the control system.

4 )Fig. 11 .

 411 The simulation of the DTP-AFPM consists of 3 stages: Stage 1: The machine operates under healthy condition; Stage 2: The phase A2 is open without reconfiguration for reference current; Stage 3: The phase A2 is open and the reference current is calculated using FTC strategy. In above descripted 3 stages, the reference torque is set to 8 Nm and the operating speed is 120 RPM. The simulated electromagnetic torque and phase current waveforms in one electrical cycle of both shift model and no shift model are given in Fig. 10 and Fig. 11, respectively. Simulation results of the no shift model: (a) Electromagnetic torque. (b) Phase current.

B

  Fig. 12. The symbol

Fig. 12 .

 12 Fig. 12. Diagram of the flux density vector components.

Fig. 10 (

 10 b) and Fig. 11(b) into the finite element analysis (FEA) model of DTP-AFPM; 3) Solve the flux density distribution in ANSYS MAXWELL, and calculate the axial stress distribution using (5); 4) Superimpose the calculated axial stress values at corresponding positions on the two sides of the rotor to get the synthetic axial stress which can be presented on a 2-D circular area as shown in Fig. 13 in which the middle lines locate at the average radius of the PMs. Furtherly, the initial circular area is unfolded to two rectangular area.

Fig. 13 .

 13 Fig. 13. Unfolding procedure from circular area to rectangular area.B. Synthetic Axial Stress in the No Shift ModelFirstly, the rotor axial force waveform of the no shift model while the machine operates in different stages is shown in Fig.14in which the symbol m θ indicates the mechanical displacement angle of the rotor. It can be seen that all the axial force waveforms fluctuate around 0 N. The synthetic axial stress distributions of the no shift model at the initial moment ( m =0 θ ) are shown in Fig. 15.Among the six subgraphs in Fig.15, there is no obvious regularity can be found. The peak value of each figure is indicated and the extremum is about 0.1MPa which is most likely due to the calculation errors.

Fig. 14 .Fig. 15 .

 1415 Fig. 14. Rotor axial force waveform of the no shift model.

17

 17

  . It can be seen that the synthetic axial stress on PMs of the shift model is variable regularly along the circumferential direction. The peak value of each figure is indicated in Fig.17

  . The comparison of the peak value of the axial stress between the no shift model and shift model is given in TABLE I. It can be found that the peak value of the shift model is higher than that of the no shift model. The detailed synthetic axial stress distributions along the middle line shown in Fig.13are given in Fig.18. The positive or negative fluctuation of the distribution curves shown in Fig.18is coincide with the position of the slot opens in stator disks. In addition, it can be noticed that the fluctuation of the axial stress on the inner PMs as shown in Fig.18(b) is higher due to the larger slot opens and smaller flux area compared with the outer PMs. It can be predicted that the possibility of the PMs breaking increases greatly compared with the no shift model as the synthetic axial stress on each PM of the shift model is uneven.

Fig. 16 . 17 .

 1617 Fig. 16. Rotor axial force waveform of the shift model.

  Fig. 18. Synthetic axial stress distribution along the middle line given in Fig. 13: (a) on the outer PMs area, (b) on the inner PMs area.

TABLE I COMPARISON

 I OF THE PEAK VALUE OF THE AXIAL STRESS (UNIT: MPA)

			No shift	Shift
	Healthy	Outer PM Inner PM	0.0758 0.0663	-0.31 -0.295
	Fault	Outer PM Inner PM	-0.0973 -0.0803	-0.377 -0.239
	FTC	Outer PM Inner PM	0.0778 -0.084	-0.288 -0.349

School of Electrical Engineering, Shandong University, China

Laboratory of Electrical Engineering and Power Electronics, University of Lille, Arts et Metiers Institute of Technology, Centrale Lille, Yncrea Hauts-de-France, ULR 2697-L2EP, France