Florian Pescher
email: florian.pescher@femto-st.fr

Benoît Piranda
email: benoit.piranda@femto-st.fr

Stephane Delalande
email: stephane.delalande@mpsa.com

Julien Bourgeois
email: julien.bourgeois@femto-st.fr

Surface Approximation by molding a shape-memory polymer on a modular robot

Keywords:

The design phase of a car development is a long and tedious process that requires a lot of trial and error. In this article we introduce a new concept aiming to make this process easier and more interactive. Our solution consist in using a modular robot together with a shape memory polymer in order to create an interactive model of a car piece. We came up with an algorithm to approximate the surface of a piece using NURBS to describe the shape and the resulting molded polymer is simulated on our simulator VisibleSim, proving the accuracy of our system.

Introduction

Nowadays during the design phase of car development, a lot of time is used creating prototypes to prepare the final style of the new vehicle. Before realizing a Computer Aided Design (CAD) of the shape, physical objects, handmade with clay, are used in order to design the most interesting parts dreamed by designers. To reduce and to optimize this fastidious work, we would like to replace physical clay with a more interactive system based on programmable matter. The matter we have in mind will be a help for CAD tools, able to display on real matter an object being designed. Interactions with the real object or the CAD model will automatically be reproduced on the other one. Thus, a designer will be able to design an object by hand or with a CAD software and restart this process as many time as needed. This matter will Fig. 1 A graph of the full system. Section A is the CAD part of the system, some ideas on how to obtain NURBS from point cloud and/or 3D model can be found in section 2. Section B is the distributed reconfiguration part of the system and detailed on the right of the figure. Explanations can be found in section 3. Section C refers to choices made by the final user, as such no further explanation is required. Section D is the polymer part of the system, details on its characteristics can be found in section 2 and the simulation implementation idea can be found in section 3 be composed of an ensemble of robots able to move by themselves around the others and equipped with a processing unit allowing them to perform calculations in order to plan their movements and achieve the desired shape. However, to realize this vision only with a modular robot, a huge amount of module is going to be necessary. To simplify the concept, we had the idea to cover the modules with a fabric able to modify its shape and recover its initial form. The most appropriate material to answer our needs is a shape memory polymer (SMP). An overview of how the complete process will work can be seen in Fig. 1 and some simulation results of the different steps are shown in Fig. 2.

Context

Robots and development context

Catoms 3D proposed by Piranda and Bourgeois in [START_REF] Piranda | Geometrical study of a quasi-spherical module for building programmable matter[END_REF] are quasi-spherical robots that can be placed in a regular Face-Centered Cubic lattice (FCC). The proposed geometry allows to connect each robots to 1 to 12 neighbors and to displace a robot to one neighbor cell of the FCC lattice by rotating it around static neighbor modules. In the present paper, we propose to use 3D catoms to define a dynamical surface mold. This mold is then used to shape the polymer suface. The displacement capabilities of 3D Catom are used to update the shape of the mold depending on an external request. Lattice system gives a position of a module with a triplet of integer coordinates. For this application, we orientate the lattice in order to define vertical columns of modules, that implies to adapt the coordinate system in the FCC lattice. Each horizontal plane ( -→ x , -→ y ) is made of stagged lines of modules, aligned along the -→ x axis. In order to get the same vector to access to each neighbor cell from a module, we define a lattice coordinate system (

- → i , - → j , - → k ) L where modules placed on - → i , - → j and - → k
axes are respectively drawn in red, green and blue. With this coordinate system, we define M as the homogeneous transformation matrix to convert lattice coordinates to world coordinates:

M =     2 × r r 0 x 0 0 √ 2 × r 0 y 0 0 r 2 × r z 0 0 0 0 1     (1)
Where (x 0 , y 0 , z 0 ) is the position of the (0, 0, 0) L cell of the lattice expressed in the world coordinate system. Every 3D catom has 12 connectors (drawn in red in Figure 3 ) used to exchange messages with neighbors. Left part of Figure 3 shows the organization of 3D catoms in the FCC lattice. In this context, we orientate the lattice in order to vertically align modules. Each vertical column is placed on a bottom connector (drawn in grey), the first module drawn in green leads the treatments of the other modules of the column. The bottom connector allows to connect every green module to a common central neighbor. These connections will be used to transfer the shape of the goal surface to every module and to synchronize the distributed treatments.

The right part of Figure 3 shows the final polymer surface deformed by the 3D catoms configuration presented on the left of the figure .   3D catoms robots are not available yet but geometrical configuration and distributed program can be evaluated in VisibleSim simulator. VisibleSim [START_REF] Piranda | Visiblesim: Your simulator for programmable matter[END_REF] is able to simulate large scale configurations of tens of thousand of robots of several shapes including 3D Catoms. The distributed program that managed robots behavior is written in C++. It allows to exchange messages between modules, actuates rotations and get events of physical interactions.

Surface Description and choice of the NURBS solution

In this section we will study the advantages and drawbacks of several shape descriptions solutions and explain our choice of using NURBS taking into account the tradeoff between precision and data size needed.

The shape description problem

In our application we first need to describe the shape we want our robots to achieve in the most efficient way, in other words using the less data possible since our robots have a limited memory. In cartesian coordinates the shape description problem can be seen as a binary classification problem : We are searching a function f such that f (x, y, z) return true if the cartesian position (x, y, z) is in the desired shape and false otherwise. In our application we are limiting ourselves to surfaces, which means that at each couple (x, y) there is a single highest point z max that describes the surface. With this assumption the shape description problem can be turned into a regression problem : We are searching for a function g such that g(x, y) returns z max . We can then obtain f from g by doing the following : f (x, y, z) returns true if z ≤ g(x, y) and false otherwise. In this section we will focus on solution that solves the regression problem.

NURBS method

NURBS [START_REF] Piegl | The NURBS Book[END_REF] or Non-Uniform Rational Basis Splines are a mathematical model used to represent curves and surfaces. Nowadays they are widely used in computer-aided design (CAD). Here we will only consider the NURBS surfaces. They are defined by a set of control points and 2 knot vectors. The control points are represented with their homogeneous coordinates ( Their 3D position (x, y, z) and a weight w). The knot vectors then define how the different control points affect the final shape. The number of knots m is defined as follow : m = n + d + 1 with n the number of control points and d the degree of the NURBS. NURBS surfaces depends on two parameters u and v, as such they have two degrees n u and n v and 2 knots vector with respectively m u and m v knots. If we note the control points P i, j with respective weight w i, j the NURBS surface S(u, v) can be computed as follow :

S(u, v) = ∑ m u -n u -1 i=0 ∑ m v -n v -1 j=0 N n u i (u)N n v j (v)w i, j P i, j ∑ m u -n u -1 i=0 ∑ m v -n v -1 j=0 N n u i (u)N n v j (v)w i, j , (u, v) ∈ [0, 1] 2 (2)
In this equation the N d j (t) are computed using the Cox-De Boor formula:

     N 0 j (t) = 1 i f t j ≤ t < t j+1 0 otherwise N d j (t) = t-t j t j+d -t j N d-1 j (t) + t j+d+1 -t t j+d+1 -t j+1 N d-1 j+1 (t) (3)
In this equation the t j are knots from the knot vector. In case some knots t j are the same we suppose 0 0 = 0.

In order to solve our regression problem with this solution we need to do some changes to the method. To find the function g we need to be able to make a clear link between the parameter (u, v) of the NURBS and the (x, y) cartesian coordinates to be able to calculate g(x, y) at any (x, y) wanted.

Multivariate polynomial interpolation

In [START_REF] Saniee | A simple expression for multivariate lagrange interpolation[END_REF] the author presents a generalization of Lagrange polynomial interpolation in multivariate case. With this method in order to obtain a polynom of m variable and of degree n you need p = n+m n points. As a generalization of Lagrange polynomial interpolation this solution present the same advantages and drawbacks : The polynom will fit the given points, however you can not really achieve a precise shape without a high degree polynom, implying a lot of points to describe the shape.

Analogy with interpolation in image processing

In image processing interpolations method [START_REF] Hs Prasantha | Comparative analysis of different interpolation schemes in image processing[END_REF] are often used, in order to scale an image for example. Those methods could be adapted to our application by doing an analogy between the height in our application and the color in the image processing context. The most widely used method in image rescaling are nearest neighbor, bilinear and bicubic interpolation. While having relatively good results in an image processing context, those method requires a lot of initial data in order to correctly reconstruct the wanted shape making it less efficient in our context.

Our choice

Now that we have studied several shape description methods we compared them in Table 1. We decided to use NURBS to describe the shape in the next steps of our work since they are widely used in CAD and can describe a shape precisely with relatively few parameters. Now that we have chosen to use NURBS to describe the shape we want to create, we need to obtain those NURBS parameters. In our system there are 2 ways of obtaining those parameters: either from a 3D Model created with a CAD software or by interacting with the robots giving back an ensemble of position, in other words a point cloud. In his thesis [START_REF] Lankinen | 3d surface modeling using point clouds[END_REF] the author describe a way of obtaining NURBS from an unorganized point cloud by first transforming it into a 3D model and then transforming it into a NURBS description. Others solutions and/or simplification in our particular application might exist but will not be the focus of this article, in the following we just consider that it is possible to obtain usable NURBS parameters from either a 3D model or a point cloud.

Polymer

Shape memory materials have the ability to change their shape and recover their original shape upon application of an external stimulus. One possible way to trigger shape memory effect is to change and increase system temperature. These material are called thermos-responsive. To obtain this ability, two conditions are required. First, switch domain as reversible thermal transition is necessary for temporary shape fixation and partial recovery. This shape memory transition allows to enable chain mobility to fix temporary shape and inversely recover permanent shape. Then, a cross-linking network determines the permanent shape to prevent chain slipping. The forming stage requires heat and the stabilization stage a temperature reduction.

A common SMP presents an extent deformation up to 800%, a density between 0,9 and 1.1 g.cm-3 and a required stress to be deformed around 1-3 MPa. Nevertheless, to limit the strength necessary to get the deformation by mems, the thickness of the foil will have to be limited consequently, a large number of transition could initiate cracks inside the material. For this reason a new material with healing properties is currently being developped and will eventually become the one used by our system.

Contribution

We are now going to detail some of the steps of our system. Left part of Fig. 1 show the complete system divided into 4 fields and the right part divide the reconfiguration process into several subfunctions. In this section we detail the distributed reconfiguration subfunctions and the polymer section of the complete system.

Reconfiguration: Flooding and robot motion functions

Transmitting the NURBS parameters, obtaining height h in (x, y) and getting position from robots are simple functions and work in a similar way. You start from the ground of pseudo catom and send messages up with either the NURBS parameters or a request for height h. Once you reach the top you send back down a message with either a "NURBS fully transmitted" response or the height h. For the getting position function you can then artificially repopulate the point cloud on the computer by adding a point every (x, y, hn * d) with d the diameter of a catom while hn * d is above ground height.

The move robots and ask for a robot in (x, y, h + d) functions are used to get from one shape to another according to the shape description we have given to the robots.

With the current work we are able to know how many robots are missing or too much in a vertical column compared to the goal shape thus creating some sinks and source of robots that we will use in future works to plan for the most efficient way to move the robots in order to achieve the goal shape.

3.2 Reconfiguration: Find wanted z max in (x, y) from NURBS

We chose to use the NURBS description method, however the usual NURBS calculation of a NURBS is a parametric surface S(u, v) = (x s , y s , z s ) with (u, v) ∈ [0, 1] 2 which makes it difficult to find our function. Since we limit our case to surfaces we can use a dichotomy method to find an approximation of g(x, y). The idea is to separate the surface into four quadrants, find the quadrant in which the couple (x, y) where we want to evaluate g(x, y) is situated and take this quadrant as our new surface. Finally we have x s almost equal to x and y s almost equal to y so we can approximate g(x, y) = z s . The algorithm is described in Algorithm 1

Algorithm 1 NURBS Dichotomy Algorithm

u 0 = 0; u 1 = 1; v 0 = 0; v 1 = 1; ε = 0.01; d = +in f inity d max =? while (d > d max ) do u = (u 0 + u 1 )/2; v = (v 0 + v 1 )/2; (x 1 , y 1 , z 1 ) = S(u -ε, v -ε); (x 2 , y 2 , z 2 ) = S(u + ε, v -ε); (x 3 , y 3 , z 3 ) = S(u -ε, v + ε); (x 4 , y 4 , z 4 ) = S(u + ε, v + ε); for i = 1 to 4 do d = (x -x i ) 2 + (y -y i ) 2 if (d < d 2 ) then quadrant = i d = √ d end if end for switch (quadrant) case 1: u 1 = u; v 1 = v; z = z 1 ; case 2: u 0 = u; v 1 = v; z = z 2 ; case 3: u 1 = u; v 0 = v; z = z 3 ; case 4: u 0 = u; v 0 = v; z = z 4 ;
end switch end while return z;

Polymer Simulation

In order to simulate the polymer we slice the polymer into an ensemble of small cubes of side length dl. We consider that each of those small cubes is submitted to the same forces than the macroscopic object itself. The full polymer is submitted to the following forces:

• Its weight : F = -m * g with g the constant of gravity and m the mass of the object • The elastic deformation defined as F dl 2 = l-dl dl * E with E Young's modulus and l the length of the contracted/extended side.

• An action of contact stopping the fall of the polymer when it hits a robot or the ground.

Since we are not considering the full polymer but an infinitesimal volume of it we can consider the volumic forces which leads to the following assumptions:

• The volumic expression of gravity is dF = -ρ * g with ρ the volumic mass • The volumic expression of elastic deformation is dF = l-dl dl 2 * E We then apply Newton's law of motion on those small cubes in order to simulate the movement of the polymer falling down on our robots in order to follow their shape.

We can see in Fig. 4 that the NURBS mathematical model seems correctly approximated by the catoms and polymer, except on the huge altitude difference on one of the sides of the object. This error can be explained by the stiffness of the polymer. We also need to take into account that in our simulation we only applied gravity to the polymer, if we had another force to attract the polymer closer to the catoms on the sides this error can be reduced experimentally. However a simple plot analysis is not enough to judge the accuracy of our method correctly so we decided to perform a statistical analysis on the data of the difference between the mathematical NURBS and the polymer. The results of this analysis can be seen in Table 2. Considering that in this experiment the diameter of a catom is 10, Table 2 shows that we can approximate really well the mathematical model on some points (really low minimum error), that the maximal error is of around 4 catoms on the side of the structure. The study of the quartiles shows that 25% of the catoms approximate the shape with an error inferior to half a catom, 50% of catoms approximate it with an error around a catom and that 75% of them shows an error of around two catoms. This analysis show an error of around 1 catom altitude in a structure contained in a 8*12*14 catoms block meaning an altitude error inferior to 10% . That leads us to wonder if this error will grow or stay constant as the structure become larger (equivalent to the catoms becoming smaller).

Scale influence

We now study the influence of the size of a catom relatively to the size of the desired shape. In order to do so we consider the previous experiment of the rear-view mirror in the 8*12*14 block as the scale 1 model and we are going to perform the same statistical analysis on different scale to see the influence of the catom size.

Table 3 shows that when the size of a catom became smaller compared to the size of the desired shape, even though the maximal error increase the values of the different quartiles remain approximately the same meaning that the relative error decrease. We can say that the smaller the catoms, the more precise our solution is.

Conclusion

In this article we introduced our new system to help the design of pieces in the car industry using a combination of modular robots and shape memory polymer. We mainly focused on finding a way to describe the desired shape introducing the NURBS dichotomy algorithm, and evaluating the accuracy of this method by simulating a polymer being molded over the catoms organized in the shape we want to design. This method proved itself to be accurate and this accuracy is improved as the catoms become smaller. Future works will focus on creating the full system and more precisely on the way to plan the catoms movement to get into the desired shape.

Fig. 2

 2 Fig. 2 The different steps of the design of a piece using our system. a) The 3D model of a rear-view mirror. b) The NURBS mathematical model representing our object. c) An ensemble of catoms taking the shape of the mirror in our simulator VisibleSim. d) The simulation of the polymer covering the ensemble of catoms.

Fig. 3

 3 Fig. 3 Left: Organization of 3D catoms along a FCC lattice. Right:Polymer surface placed over 3D catoms configuration.

Table 1

 1 Comparison between the different shape description methods

	Method	Data needed	Precision obtained	Advantages	Drawbacks
	NURBS	Polynomial order in 2 directions, 2 knot vectors, and a few control points	Good with only a few control points	Widely used in CAD. needed Good ratio precision/data	Parametric surface need to adapt to use on robots
	Multivariate Polynomial Interpolation	For a polynom of degree n in both direction need for n+2 n points	Depends on degree n. Fit the given points but bad generalization	Simple calculation of position	Need for a high degree polynom to be precise
	Nearest Neighbor	Points	Depends on number of points. Apparition of a "stair effect"	Really simple to use Need a lot of points to be precise
	Bilinear Interpolation	Points organized in a regular grid	Depends on the grid. Smoother than nearest neighbor.	Smoother than nearest neighbor	Need a lot of points Regular grid on data is a big constraint
	Bicubic Interpolation	Points organized in a regular grid	Depends on the grid. Smoother than bilinear interpolation	Smoother than bilinear interpolation	Need a lot of points Regular grid on data is a big constraint
	2.3 From 3D Model/Point Cloud to NURBS parameters

Table 2

 2 Shape accuracy statistical analysis

	Volumic mass Young Modulus min	max q1	median q3
	900 kg/m	1MPa	0.024251 43.611 5.2080 9.8923 17.232
		2MPa	0.0033179 44.741 5.0415 11.297 19.739
		3MPa	0.034316 46.927 5.2032 11.794 21.235
	1000 kg/m	1MPa	0.012972 43.213 5.1962 9.7707 16.429
		2MPa	0.044104 44.675 5.3455 10.036 18.459
		3MPa	0.011417 45.706 4.8856 11.243 20.738
	1100 kg/m	1MPa	0.014251 43.379 5.2572 10.104 17.087
		2MPa	0.0059585 44.125 5.5448 10.677 19.102
		3MPa	0.034602 46.880 5.2008 10.698 19.913

Table 3

 3 Scale influence statistical analysis

	Scale min	max q1	median q3
	1	0.014251 43.379 5.2572 10.104 17.087
	1.25 0.069377 51.736 5.6728 11.502 21.345
	1.5 0.014403 52.401 5.1813 11.364 18.264
	1.75 0.0057074 56.161 5.6487 11.455 20.554
	2	0.0061206 62.258 5.2643 11.585 18.929

Experiments

In this section we are going to test the accuracy of our system as well as the influence of the size of a catom. In order to test our system we are going to simulate the design of a car rear-view mirror with our simulator VisibleSim. We are going to run our NURBS dichotomy algorithm on the simulated catoms (see Fig. 2 c) ) and then simulate the polymer molding those catoms (see Fig. 2 d) ). The resulting polymer will then be exported into octave in order to be compared to the mathematical NURBS model. We chose to use a car rear-view mirror as a test object since this object is directly linked to the car industry and can be described as a single z = f (x, y) function. Plus the required polymer deformation to mold this shape is of only 40% so the polymer would not be torn apart during a physical test with this same shape.

Shape accuracy test

In order to test the shape accuracy of our system we plotted 3 graphs in octave : The exported data of the simulated polymer, the mathematical NURBS model and the difference between those 2 graphs.


Berk Celik
email: berk.celik@utbm.fr

Robin Roche
email: robin.roche@utbm.fr

David Bouquain

Abdellatif Miraoui

Coordinated Home Energy Management in Community Microgrids with Energy Sharing Among Smart Homes

Keywords: Energy management, smart homes, neighborhood coordination, multiagent systems, electric vehicles

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

INTRODUCTION

By enabling bi-directional data flow in the power system through advanced metering infrastructure (AMI), the smart grid enables using complex control methodologies for more efficient, economical and reliable system management [START_REF] Bozchalui | Optimal Operation of Residential Energy Hubs in Smart Grids[END_REF]. The smart grid offers all users the ability to participate activity, including for the residential sector. In this respect, classical houses become smart homes, so customers can modify their electricity consumption patterns to increase their social benefitmostly to reduce their electricity bills-through monitoring, communication and control capabilities. These features present an opportunity to customers to value their participation in a local electricity market, in interaction with the utility, and for the benefit of both sides [START_REF] Yi | Real-time opportunistic scheduling for residential demand response[END_REF]. For instance, while users can reduce their electricity bills by coordinating their actions, as will be shown in this paper, the utility can also decrease the aggregated peak load for secure and economic management of the power system.

In order to control distributed resources in the power grid, detailed data must be gathered from all parts of the power network. However, processing such large amounts of data will cause a heavy communication and computation burden for the operator [START_REF] Celik | Increasing local renewable energy use in smart neighborhoods through coordinated trading[END_REF]. Especially, the rapid penetration of renewable energy sources (RES) and electrical vehicles (EV) due to environmental awareness creates a more complex infrastructure that needs to be controlled adequately. Moreover, customers would not be pleased with sharing de-tailed information about their consumption habits and fully releasing the control of their electricity resources to another entity. Therefore, a decentralized approach where customers are the controllers of their own resources proposes a more practical, secure and privacyprotecting solution.

However, uncoordinated control can paradoxically reduce the performance of the algorithms by causing undesired issues (such as: rebound peaks, overloading, contingencies) [START_REF] Li | Automated residential demand response: Algorithmic Implications of Pricing Models[END_REF]Therefore, establishing a coordination mechanism among smart homes where the action of a user influences the decision-making of other users is necessary for achieving efficient and reliable energy management. Various studies focus on decentralized coordination among multiple smart homes. For example, in [START_REF] Guo | Decentralized Coordination of Energy Utilization for Residential Households in the Smart Grid[END_REF], an energy management algorithm is presented for an area with a load-serving-entity and multiple households with RES, storage, and controllable and non-controllable loads, in order to reduce the electricity cost of the area. In [START_REF] Mohsenian-Rad | Autonomous Demand-Side Management Based on Game-Theoretic Energy Consumption Scheduling for the Future Smart Grid[END_REF], a scheduling game is presented for consumption to reduce the electricity cost and peakto-average ratio (PAR) of the residential area. In [START_REF] Deng | Residential Energy Consumption Scheduling: A Coupled-Constraint Game Approach[END_REF], another game-theoretic approach is used to reduce the PAR of the area using appliance and EVs scheduling.

In this paper, a coordination mechanism for energy sharing among smart homes is therefore presented for the day-ahead energy management in neighborhood areas. It uses multi-agent systems (MAS), where are a suitable concept for decentralized management and enable dynamic interactions among entities. It is assumed that only users own PV and battery systems (not the utility) and they can share their residual generation as soon as it is generated (without battery) or at a later time (with battery). Furthermore, EVs are considered to provide energy for self-consumption of the smart home without sharing with neighbors. Both of the battery systems (home and vehicle) have the same charging principle where they can charge from self-generation and/or the main grid. For coordinated control, two types of entities are designed as home agents and aggregator agent. Lastly, in this study, dynamic price model is used for billing users, where the unit price is associated to the energy drawn from the main grid for the aggregated consumption of the neighborhood area. Note that grid constraints (e.g., line and transformer capacities) are not taken into account in this study.

The remainder of this paper is organized as follows. In Section 2, the system model and dynamic pricing mechanism are described. In Section 3, the problem formulation is presented. In Section 4, the proposed coordination mechanism is described. Simulation results are given in Section 5. Finally, in Section 6, the paper is concluded.

SYSTEM MODEL

The studied neighborhood model consists of U users, and each user u is connected to an aggregator through a communication link. Smart homes are connected directly to the main grid where the aggregator has no authority to control the power system, but acts as an advisor in the neighborhood to help decision-making in the smart homes. Lastly, smart homes have two-way communication ability only with the aggregator, not with each other, due to privacy concerns of the end-user.

SMART HOME ENERGY SYSTEM MODELING

In the considered smart home model, all appliances are connected to a controller (the home agent) through smart plugs where the home agent can sense and control the home appliances. However, based on user preferences and controllability, home appliances are divided into three groups: non-controllable, controllableshiftable, and controllable-interruptible appliances. Non-controllable appliances must be used whenever they are turned on, and are not controlled by the home agent. On the other hand, the consumption profile of the controllable appliances can be altered based on the user-defined and appliance operation constraints.

In the smart homes, in total, 14 types of smart appliances are modeled, where 10 are non-controllable (TV, lights, etc.), 3 are controllable-shiftable (washing machine, clothes dryer, dish washer) and one is controllable-interruptible (EV). The set of appliances is defined as L u = {1, 2, l, ..., L u } where L u is the total appliance number of home u, which can be different at each home. For example, while some users have one TV, others can have more than one. During the simulation time set T = {1, 2, t, ..., T }, the consumption profile of a smart home P c u (t) is determined as:

P c u (t) = Lu l=1 P l u (t), ∀t ∈ T (1)
where P l u (t) is the consumption profile of an appliance that consumes constant power between the start and end times t ∈ [t s r , t e r ], and nothing when t / ∈ [t s r , t e r ].

Regarding generation, residential PV systems are considered in the smart homes. However, based on the user economic situation and physical constraints of the building, users may or may not have a PV system installed, and the installed capacity is different for each home. Hence homes generation profiles are different from home to home. The output of a PV system P g u (t) is calculated by:

P g u (t) = N p u • N s u • P pv u • (G(t)/G ST C ) (2)
where N p u and N s u are the number of parallel and seriesconnected modules of the PV array, and P pv u is the rated power of the PV module. It is assumed that all smart homes are located in same geographic area, thus all the PV systems in the neighborhood area receive the same irradiance G(t) (irradiance variations are considered negligible). G ST C (t) is the irradiance value (1000 W/m 2 ) in standard test conditions (1000 W/m 2 , 25 • C).

For the energy storage system, batteries are installed only in smart homes with PV. Generally, batteries are charged when there is surplus generation, and discharged when consumption is higher than generation. However, in this study, we assume that a home control system is able to charge from the main grid, shift discharging operations to high price hours, and discharge for neighbors' consumption. Based on the determined battery injected power (given in Section 3), the battery power P b u (t) and state-of-charge (SOC) SOC(t) are determined as:

P d u /η d u ≤ P i u (t) ≤ Ṕ c u • η c u ( 3
)
P b u (t) = P i u (t) • η c u : P i u (t) > 0 P i u (t)/η d u : P i u (t) ≤ 0 (4) SOC u (t) = SOC u (t -1) + P b u (t) • △t /E b u (5) SOC min u (t) ≤ SOC(t) ≤ SOC max u (t) (6)
where η d u and η c u are the charging/discharging efficiencies, P d u and Ṕ c u are the maximum discharging/charging injection powers and SOC min u (t) and SOC max u (t) are the maximum and minimum SOC levels of the battery, and E b u is the battery capacity. Lastly, △t is the time interval between two time steps.

Although EVs are considered as an appliance in this section, the power profile of EVs shows the same characteristic as home batteries. According to same principle, the power P v,i u (t) and SOC SOC(t) v u profiles of EVs are modeled with: u is the EV battery capacity.

P v,d u /η v,d u ≤ P v,i u (t) ≤ Ṕ v,c u • η v,c u (7)
P v,b u (t) = P v,i u (t) • η v,c u : P v,i u (t) > 0 P v,i u (t)/η v,d u : P v,i u (t) ≤ 0 (8)
SOC(t) v u = SOC v u (t-1)+ P v,b u (t) • △t /E v,b u (9) SOC v,min u (t) ≤ SOC(t) ≤ SOC v,max u (t) (10

PRICE MODELING

The electric cost of customers in the neighborhood is determined using a form of dynamic pricing where the price is related to the aggregated power provided by the main grid. In this study, the same electricity price is also used for the reverse power flow from smart homes to the grid. The neighborhood electricity price is determined in two parts: a dynamic part, and a combined part. The dynamic part is modeled using a quadratic function as:

ρ(t, P n (t)) = a(t)|P n (t)| 2 + b(t)|P n (t)| + c(t) (11)
where a(t) > 0, b(t) ≥ 0 and c(t) ≥ 0 are parameters of the quadratic function, P n (t) is the aggregated net consumption of the neighborhood and ρ(t, P n (t)) is the dynamic part of the neighborhood electricity price.

After that, the dynamic part is combined with a fixedtariff d(t) which represents the wholesale market price at the upper level:

λ(t, P n (t)) = d(t) + ρ(t, P n (t)) : P n (t) > 0 d(t) -ρ(t, P n (t)) : P n (t) ≤ 0 (12
) where λ(t, P n (t)) is the electricity price scheme is used for billing users in the neighborhood. According to (12), when there is surplus generation which causes reverse flow from the neighborhood to the main grid, the electricity price will be lower than d(t), which will increase the interest of consumption at these times. Thus, the same pricing can be used for billing both types of users (consumers and producers) at the same time when reverse power flows exist.

PROBLEM FORMULATION

Home agents aim to minimize the electricity bills of their users by scheduling their controllable appliances and controlling charging/discharging operations of the home battery and EV battery. This section describes the optimization problem which is solved by the home agents.

For controlling shiftable-appliances, a scheduling interval [ ts r , te r ] is defined by the user for each appliance. Inside this interval, the home agent chooses the best time to run an appliance without jeopardizing user comfort. It is assumed that when a shiftable appliance starts operating, it cannot be stopped by the home agent until the end of its cycle. The constraint formulation for shiftable appliances is given by:

[t s r , t e r ] ∈ [ ts r , te r ] (13)
The operation of some appliances can depend on others, such as washing machines and clothes dryers. Users, logically, prefer to use a clothes dryer the after washing machine has finished its work. Therefore, this constraint is formulated by:

t s wm < t s cd -(t e wm -t s wm ) (14)
ts wm < ts cd -(t e wm -t s wm ) (15)
where wm and cd are used for indexing the washing machine and the clothes dryer.

To model the control of charging/discharging actions of the home battery system, home agents determine P b u (t) for each time interval. However, before giving the formulation, we need to mention several important assumptions. Firstly, we considered that a home battery system is able to be charged by self-generation, neighborhood generation, and the main grid. However, it cannot be charged to sell energy to the main grid while saved energy can be used for self-consumption and neighborhood consumption. Secondly, based on the modeling of the electricity profiles, the time resolution can be chosen equal to very short values (such as 1-min. resolution) for more detailed simulation. In such a model, the controller has to use a high number of inputs (such as 1440 for 1-min.) to determine the battery output at each time step, which is computationally expensive. Therefore, in this study, we define a battery control interval Z which has a lower resolution than the actual profile to reduce the battery input number in the optimization problem from T to T /Z, without changing the actual profile resolution. Based on that, P b u (t) is determined with logical inputs γ u ∈ {0, 1, 2} by:

Ru(t) = Pa(t) -P n u (t) (16)
where C total (k) is the total neighborhood cost and k is the iteration index.

Lastly, the aggregator agent determines the sold battery powers in real-time (t -domain) for the users, according to the final decision of the home agents when system reaches convergence. This results from simultaneous optimization and communication on the l -domain, as there is a possibility of mismatch existence in the t -domain with the occurrence of the condition P a (t) < P s (t). To eliminate these mismatches, the aggregator agent appliesthe proportionalsource-matching method described in [START_REF] Celik | Coordinated Energy Management Using Agents in Neighborhood Areas with RES and Storage[END_REF]:

P s u (t) = P a (t) • P s,d u (t)/P d s (t) (26)
where P s,d u (t) is the battery sold power and P d s (t) is the aggregated battery sold power at the final iteration of the decision-making process. In (26), the battery sold power of the smart homes are determined based on the ratio between P d s (t) and P s,d u (t).

RESULTS

SYSTEM SETUP

In this section, performance results of the studied case are given and compared with a baseline case where users are modeled as classic passive consumers with no communication, no coordination and no energy sharing abilities. Therefore, home batteries charge when generation is higher than consumption with self-generation, and discharge when consumption is higher than generation for self-consumption. Also, EV batteries cannot provide energy for self-consumption in the smart home. However, we assume that EV battery charging stops when the minimum required SOC is reached for next day travel in the baseline case, for a fair comparison with coordinated control.

NUMERICAL RESULTS

The daily electricity bills of users are given for the baseline and the coordinated control scenarios in Fig. 1.

With the proposed coordination mechanism, all users in the neighborhood area can reduce their electricity bills compared to the baseline scenario, although they have different types of equipments. It is a vital outcome of the coordination mechanism, because if some users cannot earn some benefit in exchange for their effort and participation, they would lose interest in controlling their resources and would turn the controller off. Hence, this situation may lead to consumer disengagement.

In Fig. 2, the neighborhood power profile from the perspective of the main grid is given. From the results, firstly, it can be seen that coordinated control achieves decreasing the aggregated peak demand of the neighborhood area by shifting the consumption of the controllable appliances to low price hours and discharging the batteries during high price hours. Especially, the effect of energy sharing can be observed around 20:00. While the home battery system is used to discharge energy for self-consumption before 20:00 in the baseline scenario, home agents kept the stored energy and discharge for self-consumption and share energy with neighbors to reduce the neighborhood consumption under the purpose of reducing the area price. Secondly, Fig. 2. Neighborhood Electricity Profile (Provided and sold energy from/to main grid). some smart homes with batteries, rather than charging with self-generation, fully charge with the aggregated surplus generation of the neighbors around 11:00, while this energy is fed back in the baseline scenario. Therefore, locally generated energy is kept inside the neighborhood area and utilized more efficiently due to the sharing capability.

Lastly, numerical results for consumption costs and peak demand powers of the neighborhood are given for both algorithms in Fig. 3. According to the total results for the neighborhood, the proposed coordination mechanism achieves 14.36% cost and 17.55% peak reduction compared to the baseline scenario.

NEXT STEPS

Although the presented coordination mechanism achieves reducing the peak consumption of the neighborhood area from the main grid, distribution system constraints are not considered in the proposed coordination and problem formulation. The effect of central generation resources (wind turbine, central battery, etc.) on the decision-making of the home agents will be investigated in future works.

CONCLUSION

In this paper, we proposed a coordination mechanism with a decentralized approach, where home agents are the decision-makers and the aggregator is the advisor in the neighborhood area. The main idea of the presented algorithm is to reduce the aggregated peak demand power of the neighborhood in addition to reducing the daily electricity bill of the users by scheduling household appliances and controlling battery (both home and EV) charging/discharging operations through a form of dynamic pricing. Agent-based modeling is used to design home and aggregator agents. Simulations results showed that the proposed coordination mechanism achieves reducing the total electricity cost and the aggregated peak consumption of the neighborhood. Moreover, all types of home users benefit from participating in the coordination mechanism, hence their receive a return for their efforts.

  ) where η v,d u and η v,c u are the EV charging/discharging efficiencies, P v,d u and Ṕ v,c u are the EV maximum discharging/charging injection powers, SOC v,min u (t) and SOC v,max u (t) are the maximum/minimum SOC levels of the EV battery, and E v,b

For

  the simulation setup, the studied neighborhood area consists of U = 20 users, where two have a home battery, PV and an EV, one has a home battery and PV, three have PV and an EV, four have just PV and one has just an EV. Electricity profiles are modeled with a 1-minute time resolution (T = 1440, △t = 1/60). Battery control and communication intervals are chosen equal to 60 minutes Z = L = 60, and price coefficients are assumed constant and taken as a(t) = 5 × 10 -5 , b(t) = 8 × 10 -4 , c(t) = 0 and d(t) = 0.16 e/kWh.Lastly, a co-simulation platform with JAVA Agent DEvelopment Framework (JADE) and MATLAB is used for agent modeling and performance evaluation. Data is exchanged between JADE and MATLAB through TCP/IP ports, by defining a unique port for each agent. Simulations are performed on a desktop computer with an Intel Core i7-3770 CPU @ 3.4 Ghz, 7.8 GB RAM and a 64-bit Ubuntu 14.04 LTS operating system.

Fig. 1 .

 1 Fig. 1. Daily Electricity Bills of the Smart Homes and Types (A = home battery + PV + EV, B = home battery + PV, C = PV + EV, D = PV, E = EV, F = None).

Fig. 3 .

 3 Fig. 3. Neighborhood aggregated (a) costs and (b) peak demand powers.

where R n u (t) is the aggregated net profile, except the net profile of user u, and P a (t) is the aggregated profile of the neighborhood. f .charge refers to full charging with P g u (t) + P o u (t) (P o u (t) is the charged power from the main grid), p.charge refers to normal charging with P g u (t), idle refers to zero power, and b.discharge refers to battery discharging with P c u (t) + R n u (t) -P g u (t). After that, the power sold by the battery discharge P s u (t) is determined with:

Based on the above formulation, the discharged battery power is used firstly for self-consumption, then sold for the neighbors consumption.

Lastly, the EV battery power is determined by using the same principle defined for the home battery system, but with two exceptions. Although EV charging is typically based on a constant current / constant voltage method, advanced methods can be implemented for EV battery control [START_REF] Paterakis | Coordinated Operation of a Neighborhood of Smart Households Comprising Electric Vehicles, Energy Storage and Distributed Generation[END_REF]. First, we assume that an EV battery is only allowed to discharge for self-consumption of a smart home (vehicle-to-home: V2H), neglecting neighborhood consumption (vehicle-to-grid: V2G). Second, there should be some energy left in the battery of the EV for the next day morning travel. According to that, EV constraints and battery power P v u (t) are determined using the same logical inputs γ v u (z) = {0, 1, 2} by:

where t v arr and t v dep are the arrival and departure times of the EV, and SOC v u is minimum required SOC for next day travel when t = t v dep . f .charge, p.charge, and idle have the same meanings and formulations than above. v.discharge refers to vehicle discharge with P c u (t) -P g u (t). The home net power profile P n u (t) is calculated by:

Finally, the objective function solved by the home agent to minimize the daily electricity bill of the user is formulated as:

s.t. eqs. (3), ( 6), ( 7), (10), ( 13), ( 14), ( 15), ( 19), (20)

It can be noted that P s u (t) is removed in (21) and added in (22), on purpose. They are needed separately during the data exchange described in Section 4 for establishing the coordination among smart homes. Note that Z is only used to ease the battery optimization problem, and not for appliances.

COORDINATION MECHANISM

In this section, the coordination mechanism is described by presenting the communication structure among neighborhood entities. Due to privacy concerns of the users, we assume that, first, home agents do not communicate with each other, and second, they use averaged data while communicating with the aggregator. We define a communication interval L where home agents can calculate the average of the actual electricity profile for each L time interval. Based on that, when home agents send a message, they convert a matrix of electricity profiles as [1 × T → 1 × T /L]. Oppositely, when home agents receive the data, they re-convert it back as [1×T /L → 1×T ] (messages are denoted using "ˆ" to represent the difference between communicated data and actual data, such that P n u (l) is the communication data of the home net electricity profile P n u (t), and l is the time index of the communication data).

At the beginning of the coordination, the aggregator agent determines the neighborhood price, the aggregated net and sold battery power profile as λ(l, Pn (l) = d(l), Pn (l) = 0 and Ps (l) = 0. With the received information, home agents minimize their objective function simultaneously, then determine and send the home net electricity profile P n u (t) and the home sold battery power P s u (t) to the aggregator. After that, the aggregator agent calculates the aggregated profile Pa (l) = U u=1 P n u (l) and the aggregated battery sold power profile Ps (l) = U u=1 P s u (l), and determines the aggregated net electricity profile with:

Next, the aggregator agent determines the neighborhood price and sends λ(l, Pn (l), Pn (l) and Ps (l) to home agents. After that, home agents run the optimization again and send the determined data back to the aggregator agent. This process continues until change in total neighborhood cost between iterations becomes negligible C total (k) -C total (k -1) ∼ = 0.

C total (k) = Pn (l) • λ(l, Pn (l)) (25)


Walid Belguerras
email: walid.belguerras@femto-st.fr

Youcef Benmessaoud
email: youcef.benmessaoud@gmail.com

Frédéric Dubas
email: fdubas@gmail.com

Kamel Boughrara
email: kamel.boughrara@g.enp.edu.dz

Mickaël Hilairet

Fréderic Dubas

Mickael Hilairet
email: mickael.hilairet@univ-fcomte.fr

F Dubas

M Hilairet

2-D Analytical Model of Conventional Switched Reluctance Machines

come

Introduction

The SRMs present many benefits for high-speed applications (e.g., electric compressor) compared with other types of machines. It can be operated at a very high-speed because it has no sliding contacts [START_REF] Gong | A novel rotor design for ultra-high speed switched reluctance machines over 1 million rpm[END_REF][START_REF] Hieu | Design of a high speed 4/2 switched reluctance motor for blender application[END_REF], and can be operate in sever conditions of temperatures. It is a competitor of a permanentmagnet synchronous machines for electric vehicle applications because of its simplicity and low cost, and its ability to operate at high-speed with low maintenance [START_REF] Takeno | Test Results and torque improvement of the 50-kW switched reluctance motor designed for hybrid electrical vehicles[END_REF][START_REF] Rain | Contributions à la commande et à la conception des machines à réluctance variable à double saillance[END_REF][START_REF] Arab | Étude et conception d'un groupe motopropulseur électrique à faibles niveau vibratoire et sonore pour véhicule électrique. Aspects 'contrôlecommande[END_REF].

In the literature, we find different methods of electromagnetic modeling of electric machines; semi-analytical modeling based on the magnetic equivalent circuit (or permeance network) [START_REF] Benlamine | Modeling of an axial-flux interior PMs machine for an automotive application using magnetic equivalent circuit[END_REF][START_REF] Benlamine | Nonlinear Adaptive Magnetic Equivalent Circuit of a Radial-Flux Interior Permanent-Magnet Machine Using Air-Gap Sliding-Line Technic[END_REF], subdomain method in linear conditions (i.e., infinite permeability of the iron parts) [START_REF] Dubas | Analytical Solution of the Magnetic Field in Permanent-Magnet Motors Taking Into Account Slotting Effect: No-load Vector Potential and Flux Density Calculation[END_REF][START_REF] Lubin | 2-D Exact Analytical Model for Surface-Mounted Permanent-Magnet Motors With Semi-Closed Slots[END_REF][START_REF] Lubin | Exact Analytical Method for Magnetic Field Computation in the Air Gap of Cylindrical Electrical Machines Considering Slotting Effects[END_REF], and the exact subdomain method taken into account the iron permeability [START_REF] Dubas | New Scientific Contribution on the 2-D Subdomain Technique in Cartesian Coordinates: Taking into Account of Iron Parts[END_REF][START_REF] Dubas | New Scientific Contribution on the 2-D Subdomain Technique in Polar Coordinates: Taking into Account of Iron Parts[END_REF][START_REF] Roubache | New Subdomain Technique for Electromagnetic Performances Calculation in Radial-Flux Electrical Machines Considering Finite Soft-Magnetic Material Permeability[END_REF][START_REF] Yahia | A 2-D Exact Subdomain Technique in Switched Reluctance Machines Taking Into Account of Finite Soft-Magnetic Material Permeability[END_REF][START_REF] Ben Yahia | Two-Dimensonal Exact Subdomain Technique of Switched Reluctance Machines with Sinusoidal Current Excitation[END_REF]. In addition, we find analytical methods based on multi-layers [START_REF] Sprangers | Semi-Analytical Framework for Synchronous Reluctance Motor Analysis Including Finite Soft-Magnetic Material Permeability[END_REF][START_REF] Sprangers | Magnetic Saturation in Semi-Analytical Harmonic Modeling for Electric Machine Analysis[END_REF][START_REF] Djelloul-Khedda | Nonlinear Analytical Prediction of Magnetic Field and Electromagnetic Performances in Switched Reluctance Machines[END_REF][START_REF] Djelloul-Khedda | Semi-Analytical Magnetic Field Predicting in Many Structures of Permanent-Magnet Synchronous Machines Considering the Iron Permeability[END_REF] or elementary subdomains for the local saturation effect [START_REF] Roubache | Elementary subdomain technique for magnetic field calculation in rotating electrical machines with local saturation effect[END_REF][START_REF] Roubache | Technique en Sous-Domaines Élementaires dans les Machines Asynchrones à Cage d'Écureuil : Saturation Magnétique Locale & Courants de Foucault dans les Barres[END_REF]. The analytical methods cited previously give an accurate electromagnetic result compared with numerical calculations, with reduced computation time. In this paper, we will present a comparison study of 2-D electromagnetic performances between the developed linear model based on the subdomain method in linear conditions. In order to analyse the validity of the developed model, the results has been compared with those computed numerically using FEMM [START_REF] Meeker | Finite Element Method Magnetics[END_REF] in linear and no-linear conditions.

Analytical model

The analytical model based on the subdomain method in linear conditions is given in [START_REF] Dubas | Analytical Solution of the Magnetic Field in Permanent-Magnet Motors Taking Into Account Slotting Effect: No-load Vector Potential and Flux Density Calculation[END_REF][START_REF] Lubin | 2-D Exact Analytical Model for Surface-Mounted Permanent-Magnet Motors With Semi-Closed Slots[END_REF][START_REF] Lubin | Exact Analytical Method for Magnetic Field Computation in the Air Gap of Cylindrical Electrical Machines Considering Slotting Effects[END_REF]. In order to simplify the model, we have considered the following assumptions: -End-effects are neglected, i.e.,

 

0; 0; z A  A ;
-Eddy-currents effects in all materials are neglected; -Current density in the stator slots has only on component along the z-axis, i.e.,   0; 0; z JJ  ; -The slots have a radial sides; -The relative permeability is considered infinite for the iron parts (i.e., the saturation effect is neglected). The schematic representation of the studied 8/6 conventional SRM is shown in Fig. 1. In developing the 2D analytical model, a magnetic vector potential formulation is used in polar coordinates. It consists of solving the partial differential equation [START_REF] Farlow | Partial Differential Equations for Scientists and Engineers[END_REF] due to Maxwell's equations in each domain of the studied machine. The equations to be solved in each region are given by

             (Region I), (1a)
             (Region j th ). ( 1b
) 0 2 2 2 11 zi zi zi zi A A A J rr rr                 (Region i th ). ( 1c
)
Solving of PDEs given by (1) allows to obtain the general solution of z A in each domain [START_REF] Dubas | Analytical Solution of the Magnetic Field in Permanent-Magnet Motors Taking Into Account Slotting Effect: No-load Vector Potential and Flux Density Calculation[END_REF][START_REF] Lubin | 2-D Exact Analytical Model for Surface-Mounted Permanent-Magnet Motors With Semi-Closed Slots[END_REF][START_REF] Lubin | Exact Analytical Method for Magnetic Field Computation in the Air Gap of Cylindrical Electrical Machines Considering Slotting Effects[END_REF]. The integration constants are determined by using a Fourier series expansion of z A in each region and the boundary conditions (BCs) [START_REF] Dubas | Analytical Solution of the Magnetic Field in Permanent-Magnet Motors Taking Into Account Slotting Effect: No-load Vector Potential and Flux Density Calculation[END_REF][START_REF] Lubin | 2-D Exact Analytical Model for Surface-Mounted Permanent-Magnet Motors With Semi-Closed Slots[END_REF][START_REF] Lubin | Exact Analytical Method for Magnetic Field Computation in the Air Gap of Cylindrical Electrical Machines Considering Slotting Effects[END_REF][START_REF] Dubas | New Scientific Contribution on the 2-D Subdomain Technique in Cartesian Coordinates: Taking into Account of Iron Parts[END_REF][START_REF] Dubas | New Scientific Contribution on the 2-D Subdomain Technique in Polar Coordinates: Taking into Account of Iron Parts[END_REF][START_REF] Roubache | New Subdomain Technique for Electromagnetic Performances Calculation in Radial-Flux Electrical Machines Considering Finite Soft-Magnetic Material Permeability[END_REF][START_REF] Yahia | A 2-D Exact Subdomain Technique in Switched Reluctance Machines Taking Into Account of Finite Soft-Magnetic Material Permeability[END_REF][START_REF] Ben Yahia | Two-Dimensonal Exact Subdomain Technique of Switched Reluctance Machines with Sinusoidal Current Excitation[END_REF][START_REF] Sprangers | Semi-Analytical Framework for Synchronous Reluctance Motor Analysis Including Finite Soft-Magnetic Material Permeability[END_REF][START_REF] Sprangers | Magnetic Saturation in Semi-Analytical Harmonic Modeling for Electric Machine Analysis[END_REF][START_REF] Djelloul-Khedda | Nonlinear Analytical Prediction of Magnetic Field and Electromagnetic Performances in Switched Reluctance Machines[END_REF][START_REF] Djelloul-Khedda | Semi-Analytical Magnetic Field Predicting in Many Structures of Permanent-Magnet Synchronous Machines Considering the Iron Permeability[END_REF][START_REF] Roubache | Elementary subdomain technique for magnetic field calculation in rotating electrical machines with local saturation effect[END_REF][START_REF] Roubache | Technique en Sous-Domaines Élementaires dans les Machines Asynchrones à Cage d'Écureuil : Saturation Magnétique Locale & Courants de Foucault dans les Barres[END_REF]. The linear systeme can be written as

      A X B , (2)
where: -

 

A is the square matrix of the integration constants obtained by BCs of dimension QQ  with

    4 1 1 r s Q N Q M Q K        in which s Q & r Q
represents respectively number of the stator and rotor slots, and N , M & K represents respectively the finite number of spatial harmonics terms in various regions; -  X is the vector of unknowns (integration constants to Using the geometrical and physical parameters given in Table 1, we have calculated the computing time necessary for obtaining the vector   X (viz., t=58.85 s for 31 position), and the computing time necessary for making the mesh and analyse on FEMM (viz., t=151.33 s for 31 positions. The auto-mesh is used).

Simulation results

The analytical expression of electromagnetic torque and the method for calculating the flux are given in [START_REF] Lubin | 2-D Exact Analytical Model for Surface-Mounted Permanent-Magnet Motors With Semi-Closed Slots[END_REF][START_REF] Lubin | Exact Analytical Method for Magnetic Field Computation in the Air Gap of Cylindrical Electrical Machines Considering Slotting Effects[END_REF][START_REF] Roubache | New Subdomain Technique for Electromagnetic Performances Calculation in Radial-Flux Electrical Machines Considering Finite Soft-Magnetic Material Permeability[END_REF][START_REF] Yahia | A 2-D Exact Subdomain Technique in Switched Reluctance Machines Taking Into Account of Finite Soft-Magnetic Material Permeability[END_REF][START_REF] Ben Yahia | Two-Dimensonal Exact Subdomain Technique of Switched Reluctance Machines with Sinusoidal Current Excitation[END_REF][START_REF] Sprangers | Semi-Analytical Framework for Synchronous Reluctance Motor Analysis Including Finite Soft-Magnetic Material Permeability[END_REF][START_REF] Sprangers | Magnetic Saturation in Semi-Analytical Harmonic Modeling for Electric Machine Analysis[END_REF][START_REF] Djelloul-Khedda | Nonlinear Analytical Prediction of Magnetic Field and Electromagnetic Performances in Switched Reluctance Machines[END_REF][START_REF] Djelloul-Khedda | Semi-Analytical Magnetic Field Predicting in Many Structures of Permanent-Magnet Synchronous Machines Considering the Iron Permeability[END_REF][START_REF] Roubache | Elementary subdomain technique for magnetic field calculation in rotating electrical machines with local saturation effect[END_REF][START_REF] Roubache | Technique en Sous-Domaines Élementaires dans les Machines Asynchrones à Cage d'Écureuil : Saturation Magnétique Locale & Courants de Foucault dans les Barres[END_REF]. We have used the parameters of SRM given in [START_REF] Hannoun | Etude et mise en oeuvre de lois de commande de la machine à réluctance variable à double saillance[END_REF] for our comparison study. Fig. 2 shows the equipotential lines of z A in the machine due to phase A (IA=20 A) at 180° rotor position obtained by numerical model. Fig. 3 shows the radial and tangential magnetic flux density in the middle of the air-gap by feeding only phase A (I=10 A) for rotor position 45°. Fig. 4 shows a comparison between the numerical and analytical result of flux per phase due to phase A. In linear (I=10 A), the numerical and analytical results are in well agreement. In no-linear, the relative error between analytical and numerical results is 46.56 %. The mutual flux between phase A and others phases is shown in Fig. 5. In no-linear, the analytical and numerical results present the 49.12 % relative error. The mutual flux between phases A-C is null because the opening between this phases is 𝜋/2. Fig. 6 shows the comparison between analytical and numerical results of static electromagnetic torque due to phase A. It can be seen that the analytical model gives the accurate results in linear. In no-linear, we have a 42.15 % relative error between the analytical and numerical results. This important relative error due to no considering the relative permeability of iron parts in the analytical model.

Conclusion

In this paper, we have presented an analytical model based on the subdomain method of 8/6 conventional SRM in linear conditions. The analytical results are in excellent agreement with numerical ones. However, in no-linear conditions (I=50 A), the developed model overestimates the electromagnetic performances with a maximum error of about 46 %. However, by taking into account this error, this model can be usefully used in optimization procedure where saving of computation time is required. In order to obtain more accurate performances in no-linear conditions, this model can be extended to the saturation case (taking into account the characteristic of relative permeability of iron parts). It will be the object of our future works.

Fig. 1

 1 Fig. 1 Example of 8/6 conventional SRM.

B- 11 s

 11 is the vector of electromagnetic sources terms with dimension 1 Q  . The vector   X can divided in 3 parts as following: -Part 1 is air-gap (i.e., Region I) with dimension 41 N  ; -Part 2 is the rotor slots (i.e., Region j th ) with dimension Part 3 is the stator slots (i.e., Region i th ) with dimension   QK    .

Fig. 2 Fig. 3

 23 Fig. 2 Equipotential lines of zA due to phase A for IA=20 A and 180° rotor position.

Fig. 4

 4 Fig. 4 Waveform of flux per phase due to phase A.

Fig. 5 Fig. 6

 56 Waveform of mutual flux obtained by feeding phase A in: (a) linear (IA=10 A), and (b) no-linear (IA=50 A). Waveform of the static electromagnetic torque due to phase A in: (a) linear (IA=10 A), and (b) no-linear (IA=50 A).

Fig. 7

 7 Fig.[START_REF] Benlamine | Nonlinear Adaptive Magnetic Equivalent Circuit of a Radial-Flux Interior Permanent-Magnet Machine Using Air-Gap Sliding-Line Technic[END_REF] The maximum electromagnetic torque due to phase A for different values of current.

Fig. 8

 8 Fig.[START_REF] Dubas | Analytical Solution of the Magnetic Field in Permanent-Magnet Motors Taking Into Account Slotting Effect: No-load Vector Potential and Flux Density Calculation[END_REF] The maximum self-inductance due to phase A for different values of current.

Fig. 7

 7 Fig.7shows the comparison between analytical model and (no-)linear numerical results of the maximum electromagnetic torque due to phase A for different values of current (0-100 A). The results of the max self-inductance due to phase A for different values of current (0-100 A) obtained by analytical model, and (no-)linear numerical model are shown in Fig.8. It can be seen that the analytical model gives a good results with no-linear numerical result in limited current rang (I ≤ 30 A for max torque, and I ≤ 20 A for max self-inductance).

Table 1

 1 Parameters of 8/6 conventional SRM Parameter Value Q

		opening	19.8°
	θ rso Rotor slot opening	39.3°
	θ tr Rotor teeth opening	20.6°
	N Number of harmonic in air-gap	150
	M Number of harmonic in rotor slots	150
	K Number of harmonic in stator slots	150
	P	Pitch used in calculation	12°

s Number of stator slots 8 Q r Number of rotor slots 6 R sh Radius of the shaft 11.5 mm R ry Radius of the rotor yoke 20.9 mm R r Radius of the rotor surface 33.9 mm R s The stator bore radius 43 mm R sy Radius of the stator yoke 53.4 mm R out Outer radius of the stator 71.5 mm L Axial length 125 mm θ 𝑠𝑠𝑜 Stator slot opening 25.2° θ 𝑤 Opening of a slot coil 12.6° θ ts Stator teeth


E Saint Ellier

Y Bailly

L Girardot

D Ramel

P Nika

ACOUSTIC STREAMING MEASUREMENTS BY PIV

Keywords: Main subjects: flow visualization, thermoacoustics Fluid: pulsed flow Visualization method(s): particle image velocimetry Other keywords: image processing, acoustic streaming

The acoustic streaming is a second-order steady flow, which is superimposed on the dominant acoustic velocity. It is induced by the nonlinearities of the acoustic propagation. These non linear phenomena are present in oscillating flows and generally result in a efficiency degradation of thermoacoustic devices. In this paper an exploration of the acoustic velocity field by PIV enabled to highlight and quantify these secondary flows. The measurements performed in phase synchronisation and by scanning the temporal period give us information on streaming profiles and allow plotting the acoustic velocity over time.

INTRODUCTION.

The study of acoustic resonators can lead to a variety of subjects, from music [START_REF] Bork | Practical tuning of xylophone bars and resonators[END_REF] to photoacoustic spectrometry [START_REF] Tavakoli | Design, simulation and structural optimization of a longitudinal acoustic resonator for trace gas detection using laser photoacoustic spectroscopy (LPAS)[END_REF], [START_REF] Miklós | Application of acoustic resonators in photoacoustic trace gas analysis and metrology[END_REF] through thermoacoustics [START_REF] Penelet | Comportement transitoire d'un moteur thermoacoustique annulaire[END_REF], [START_REF] Yazaki | Traveling Wave Thermoacoustic Engine in a Looped Tube[END_REF], [START_REF] Duffourd | Etudes expérimentales d'un stack de réfrigérateur thermoacoustique : mesure de vitesse acoustique par PIV et évolution temporelle du gradient de température[END_REF]. The present study is motivated by this last discipline. Thermoacoustics is actually a relatively new research field in physics which started in 1980s at the Los Alamos National Laboratory (LANL) in the USA. It is based on the reverse conversion between acoustic and thermal energy. The systems use an environmentally friendly working medium (noble gas) in a Stirling-like cycle, and contain no moving parts which give them a high reliability as well as a long life span. It is a clear, reliable energy and further development of these systems will lead to attractive renewable energy options. The most basic systems consist of an acoustic resonator within which is located a stack equipped with two heat exchangers at its ends. In the case of an engine, if the temperature difference maintained by the exchangers between the stack sides is sufficient, it can create an acoustic wave. In the case of a receiver, it is the fluid oscillations created by an acoustic generator that will allow pumping heat along the stack from the cold area to hot area. Most applications involve heating cooling or electricity production. To be more specific thermoacoustic systems can for example be applied as heat pumps for domestic applications or industrial applications like liquefaction of natural gas or solar driven cooling systems. If this new technology, clean and reliable is a clear interest, it still remains a challenge. Indeed in thermoacoustic engines and heat pumps significant thermal losses are present due to non-linear effects. An important class of these non-linear effects is streaming flows, which are secondorder steady flows, superimposed on the first-order acoustic oscillations. Several types of streaming exist: Gedeon streaming, Rayleigh streaming or jet streaming. One type of the acoustic streaming which is always associated with a standing-wave resonator is called Rayleigh streaming. It refers to the streaming motion that is induced by acoustic standing waves in thermoviscous fluids between plane walls. It actually originates from the action of Reynolds stresses developing in the Stokes boundary layer (boundary layer near a solid wall in a viscous and oscillating fluid) [START_REF] Riley | Acoustic streaming[END_REF]. These flows create steady mass-fluxes that transport heat away from the hot heat exchanger to the cold heat exchanger of a thermoacoustic engine, without contributing to energy conversion, or in case of a refrigerator in the opposite direction of the heat pumping, leading in both cases to a degradation of performance. The understanding of these non-linear effects will therefore be important in improving system performances. It is thus necessary to develop a methodology to provide insight into these phenomena. Although acoustic streaming has been widely investigated, few works on experimental measurements were performed. Studies generally confine to local pressure and acoustic velocity measurements. Evaluating experimentally a non-linear effect of the acoustic wave such as acoustic streaming is a real challenge in metrology and signal processing. In the field of thermoacoustics Poignand [START_REF] Poignand | Réfrigérateur thermoacoustique : etude du système compact et du comportement transitoire[END_REF] used LDA to measure the particle velocity inside a thermoacoustic refrigerator composed of one cubic cavity with two or four loudspeakers attached to the walls. Laser Doppler Anemometry has a very good time resolution and therefore gives very good results when one wishes to have access to the local temporal evolution of the velocity. In the previous study PIV was used jointly with LDA. Unlike LDA, the Particle Image Velocimetry enables instantaneous measurements of the spatial field. For that matter PIV was also used by Debesse and al. [START_REF] Debesse | Adaptation de la vélocimétrie par images de particules à l'analyse des effets non linéaires en thermoacoustique[END_REF] to measure acoustic streaming in thermoacoustic systems. The singular value decomposition (SVD) allows them to reconstruct the acoustic signal despite the sub sampling without requiring a phase reference. The mean-field study, conducted on an almost 7m-long resonator, shows streaming flows with velocity amplitudes between 1 and 8 cm/s. Yet the authors remain reticent as for the repeatability of these second-order flows. Studies were also carried out on simple acoustic resonators. Thompson and Atchley [START_REF] Thompson | Simultaneous measurements of acoustic and streaming velocities in a standing wave using laser Doppler anemometry[END_REF] used LDA to study the acoustic velocity as well as the streaming velocity inside a cylindrical standing-wave resonator filled with air. They showed that the magnitude of the streaming is influenced by the dependence of fluid viscosity with temperature. This mean of measurement was also used by Sourice and al. [START_REF] Sourice | Mesures de vitesses acoustiques en temps-réel par LDV[END_REF]. In fact they developed an algorithm for signal processing allowing real-time estimation of the acoustic velocity. Eventually this algorithm was implemented on a LDA system. Their method was validated on a simple case, in a Plexiglas tube excited by an electrodynamic loudspeaker. Gazengel and Poggi [START_REF] Gazengel | Measurement of acoustic particle velocities in enclosed sound field: Assessment of two Laser Doppler Velocimetry measuring systems[END_REF] evaluated two LDA measurement systems. Reference velocities were defined from pressure measurements combined with a model of sound propagation. The equipment was then put in place to minimize errors on the velocity reference. These velocities are used subsequently to validate the results of the LDA. PIV and LDA are not the only means of measurements used to measure acoustic velocities. Huelsz and Lopez-Alquicira [START_REF] Huelsz | Hot-wire anemometry in acoustic waves[END_REF] developed a method using a hot wire probe suitable for measuring sound waves with low acoustic amplitude. Their measurements within a quarter-wave resonator nevertheless outline the presence of streaming. The process of hot-wire probe has the advantage of having relatively good spatial and temporal resolution. Also the implementation is generally simple. Besides compared with the PIV method, the measurement requires no seeding and therefore allows quick installation at low cost. However this method requires introducing a probe into the flow, which remains essentially intrusive and sensitive to changes of temperatures that can potentially twist the result. The measure is also local which requires to cover a large number of points to get a representation of the velocity field. Few two-dimensional experimental investigations on the non-linear phenomena of acoustic streaming inside a standing wave resonator have been done. Nabavi and al. [START_REF] Nabavi | Measurement of the acoustic velocity field of non linear standing waves using the synchronized PIV technique[END_REF] applied the PIV to the standing acoustic waves developing in a 7cm² -square section resonator with air at atmospheric pressure. Measurements of the acoustic field were first performed for quasi-linear conditions in order to validate them by comparison with the results of an analytical model. The same authors [15] conducted experimental measurements inside a Plexiglas channel of 4cm² -square cross-section filled with air (the setup is similar to what was used in [START_REF] Nabavi | Measurement of the acoustic velocity field of non linear standing waves using the synchronized PIV technique[END_REF]) to analyze the formation of regular and irregular acoustic streaming patterns. When the excitation amplitude of the acoustic driver increases beyond a certain limit, regular streaming flow patterns are distorted to an irregular flow structure.

EXPERIMENTAL SETUP.

Thermoacoustic bench.

In this paper the system consists of a resonator made of stainless steel closed at one of its ends and equipped with an acoustic driver to the other. The inside diameter of the resonator is 56.3mm. The resonator is made of six segments, in which is included the PIV module prepared for flow illumination with laser light and imaging. It has three windows of visualization to let the laser beam go through and a prism to deflect the laser sheet onto the upper window into the resonator (see Fig. 1). The cell was manufactured with square section and smooth angles. At both ends it allows a smooth section change from square to circle to fit the resonator, without modifying the cross section area. The optical elements for visualization have been treated against-reflection for the wave length 532nm (see the pink color on Fig. 1), which is the laser wavelength. The shaker, acting as acoustic source, (LDS -Ling Dynamic Systems-V450/1 -PA 500L) was chosen for reliability reasons considering its working conditions and thus matching those resources to provide as much as possible for monochromatic waves without getting measurements altered by parasitical thermal effects. It provides a mechanical June 25-28, 2012, Minsk, Belarus ISFV15 -Minsk / Belarus -2012 vibration which amplitude can be selected by modifying the piston stroke. Consequently the piston displacement and the operating frequency can be adjusted. The piston diameter matches the resonators inside diameter of 56.3mm. Its stroke varies from 5 to 18mm. The image acquisition is then achieved with a CCD camera (TSI PIvCam 13-8) of 1024x1248 pixels. The images from the camera are recorded at a frequency of 3.63 Hz. The camera is connected to a synchronizer (LASERPULSE Synchronizer -TSI Model 610034) allowing triggering the image acquisition with the control signal of the shaker. Consequently the PIV measurements are synchronized with the acoustic wave. An aerosol generator (TOPAS ATM 210) is used to generate the seeding mist. The particles are made of DEHS (Di-2-Ethylhexyl-Sebacat). The size of the droplets produced is of the order of 0.3  m and the droplets evaporate completely after 4 hours thus preventing clogging. Several valves control the seeding and its injection in the resonator. The plenum chamber of the seeding system can withstand a pressure up to 10 bars and allows the homogenization of the gas-particles prior to injection into the measuring zone [16]. Pressure and temperature acquisition. In addition to the PIV system, the resonator is equipped with pressure and temperature sensors (see Fig. 3). The pressure transducers are Kulite transducers (WCT-312 and XTEL-190) connected to an acquisition board National Instrument (Ni SCX1-1000). A K type thermocouple provides the average temperature within the PIV cell. The thermocouple allows knowing the temperature inside the resonator, temperature from which can be calculated the acoustic wave velocity thereby deducing the resonance frequency.

The set of signals is captured at a sample rate of 15 KHz. The "drive ratio"-defined as the ratio of the maximum amplitude of acoustic pressure over the mean pressure-is calculated from pressure measurements taken at the resonator's plug (closed end).

Post processing

Image acquisition is done using Insight 3G software that controls various parameters: -LPD (ie, its Laser Pulse Delay): time between the acquisition trigger order and the first laser flash = 220s -The t= time between two laser flashes = 25s -PIV-Cam frequency: acquisition frequency of each pair of images (here 3.63 Hz, which is the maximum allowed by the camera). 50 pairs of images are recorded for each phase knowing that twenty-one phases are studied on the acoustic cycle. The trigger of PIV occurs on the rising edge of the signal generated by the synchronizer from the acoustic signal. The image quality has an influence on the correlation performed by the PIV software (Insight 3G). Appropriate treatment/processing is then operated on different images: the noise is first subtracted from each snapshot, and then several operations are performed on the instantaneous velocity vectors fields.

1. The instantaneous velocity fields are statistically averaged. A mean velocity field is then created for each phase.

2. To overcome the side effects, only a portion of the image is kept. 3. For streaming velocity fields, the fields are averaged over the acoustic signal. 4. Velocity profiles and a temporal representation of the acoustic velocity variation over a time period can be then deduced from the velocity fields.

RESULTS

Drive ratios (Dr) from 3% to 6% were obtained with the following experiment conditions:

-Resonance frequency: 25Hz -Atmospheric pressure -Fluid: air The detailed experiment conditions are summarized in the following The next figure (Fig. 4) presents the temporal variation of the axial velocity inside the resonator for four different drive ratios. The velocity, averaged spatially, is shown over time in ms, for one period. The shape of the curves is far from being sinusoidal and is actually much distorted. In addition the higher velocities, presented in the second half of the period, are revealing a strong dissymmetry of the axial velocity, which drives away the particles from the shaker. This pattern is most likely explained by the demonstration of a non linear phenomenon and probably the presence of superior harmonics, distorting the sinus. Comparatively the piston amplitude, presented on Figure 5, follows a perfect sine curve, which indicates that the presence of superior harmonics in the flow isn't related to the acoustic driver. Figure 6 presents the streaming velocity profiles inside the PIV cell, on a plane located at the center of the resonator (xy plane). In the y-axis is the height of the resonator in mm and in the x-axis is the axial velocity in m/s. The shape, well defined at high drive ratio (see Dr = 5.5%), gives a vortex-like streaming velocity or what is called the Rayleigh streaming (Rayleigh cell) with opposite velocities at the center and on the wall of the resonator, which was expected [16]. Yet the resolution of camera doesn't allow to have the complete shape of the cell velocity profile since detailed data close to the walls are not available. The amplitudes are similar to what was found in the work of P. Debesse [17]. However the 3-D characteristic of the streaming flow was not studied here and has to be taken into account if we want the fully characterize the streaming.

CONCLUSION

The work presented is part of an experimental exploration phase of acoustic streaming in thermoacoustic systems. This first phase has allowed to adapt the technique of PIV to the specific character of Rayleigh streaming considering the low intensity of this order-two phenomenon that under certain conditions is comparable to the accuracy of current equipment.

Measurements were also performed under specific conditions dictated by the imperatives of geometry: low pressure (air at atmospheric pressure) and isothermal flow. In addition, the three-dimensionality of the flow streaming is not studied here and should nevertheless be taken into account if one wishes to fully characterize the fluidic phenomena. However the results allowed revealing the nonlinear character of the acoustic wave in conditions close to the thermoacoustic (drive ratio around 5%). Thus, in addition to the higher harmonics that distort the sound signal, streaming profiles measured in the xy plane appear to correspond to representations of Rayleigh streaming observed in pipes of length suitable to the establishment of a standing wave. Currently it is planned to perform further experimental measurements on a Plexiglas channel of a square section allowing optical visualisation along the whole resonator.

Fig. 1

 1 Fig. 1 Views of the PIV cell, treated against reflections

Fig. 2

 2 Fig. 2 Schematic of experimental apparatus and instrumentation

Fig. 3

 3 Fig. 3 Acquisition process pressure and temperature signals

Fig. 4 Fig. 5

 45 Fig.[START_REF] Penelet | Comportement transitoire d'un moteur thermoacoustique annulaire[END_REF] Time variation of the axial velocity over a period for different drive ratios, ×: Dr=3.3%, : Dr=4.2%,: Dr=4.8%: Dr=5.5%

Fig. 6

 6 Fig.[START_REF] Duffourd | Etudes expérimentales d'un stack de réfrigérateur thermoacoustique : mesure de vitesse acoustique par PIV et évolution temporelle du gradient de température[END_REF] Streaming velocity profiles for different drive ratios, ×: Dr=3.3%, : Dr=4.2%,: Dr=4.8%: Dr=5.5%

table . June 25-28, 2012, Minsk, Belarus

 .

		Test1	Test2	Test3	Test4
	P m abs. (bar)	1,0146 1,0108 1,0165 1,0140
	P 1 (mbar)	34	42	49	56
	Drive ratio P 1 /P m (%) 3,3	4,2	4,8	5,5
	T m (°C)	21,7	21,9	21,5	21,8
	Piston stroke (mm)	7,07	8,48	9,90	11,3
	T (period in s)	0,040	0,040	0,040	0,040
	f (Hz)	25	25	25	25
	L (m)	6,85	6,85	6,85	6,85
	Acoustic driver	Shaker Shaker Shaker Shaker
	Tab. 1 Experiment conditions of the tests conducted in air
	With P m : mean pressure in the resonator
	P 1 : acoustic pressure
	L: resonator length
	ISFV15 -Minsk / Belarus -2012

ACKNOWLEDGMENTS:

This work is supported by the European Union under the FP7. Contract No. 226415 -www.thatea.eu

15. Nabavi M., Siddiqui K., Dargahi J., Analysis of regular and irregular acoustic streaming patterns in a rectangular