S Zafar Satie

G Camilleri
email: guy.camilleri@irit.fr

V Maurya

H Ben

Ahmed Blavette Satie

M-P Gleizes

Adaptive Multi-Agent System and Mixed Integer Linear Programming Optimization Comparison for Grid Stability and Commitment Mismatch in Smart Grids

Keywords: Multi-Agent Systems, Mixed Integer Linear Programming, Optimality Gap, Electric Vehicles, Balance Responsible Party, Distribution Network Stability

Existing electrical networks are going through a transition and distributed energy resource, if not managed properly, can hinder this transition. Uncontrolled introduction of photovoltaics and electric vehicles in distribution networks would lead to substantial issues such as commitment mismatches, line congestions, voltage deviations, etc. This paper presents the use of a classical approach, mixed integer linear programming optimization, and a novel approach, adaptive multi-agent system, to solve the highlighted distribution side challenges by utilizing electric vehicles' storage capacity. This comparison serves as a great tool to benchmark the performance of the under-development adaptive multi-agent system methodology.

I. INTRODUCTION

Electrical networks, at first designed for unidirectional power flows, are facing stability and resilience challenges due to introduction of modern distributed energy resources (DERs) such as photovoltaics (PVs) and electric vehicles (EVs). The European Union (EU) 2030 climate and energy framework [START_REF]2030 climate & energy framework[END_REF] and projections by the International Energy Agency (IEA) [2] present a future with appreciable amount of such DERs connected to the distribution network. This certainly ensures a smoother path towards carbon-free energy systems but can also cause disruptions in the existing electrical infrastructure by bringing techno-economic challenges.

On the technical side, congestions in the network, voltage deviations as well as a need for energy balancing due to imperfect forecasts can occur. The financial responsibility of balancing is managed by Balance Responsible Parties (BRPs). BRPs submit a daily balance schedule to the Transmission System Operator (TSO) and manage energy balance within their portfolio called the Balance Perimeter (BP). This schedule is defined on a sub-hourly basis (soon to be harmonized to a quarter-hourly basis in Europe), known as the imbalance settlement period. If the BRP has consumed less energy (i.e., produced more energy) than what was planned, it would be compensated by the TSO and vice versa. To mitigate the commitment mismatch cost, a BRP may use, through an aggregator, flexibles entities (e.g., EVs). However, the BRP orders must not create local issues on the network, such as congestions. Indeed, the problem resides in the interactions between different actors managing different roles at different levels and having potentially contradictive objectives. Hence, with an increasing amount of DERs, it is necessary to mitigate these challenges.

Different solutions have been proposed to ensure an efficient coordination while maintaining grid stability. Grid reinforcement solutions are somewhat limited due to a significant infrastructure cost, therefore making flexibility solutions more promising. One such approach is the application of Multi-Agent System (MAS) concepts to manage networks. Hu et al. [START_REF] Hu | Multiagent-based modelling for electric vehicle integration in a distribution network operation[END_REF] have used swarm intelligence to model a MAS for controlling a distribution network with EVs. Lauri et al. [START_REF] Lauri | Managing Power Flows in Microgrids Using Multi-Agent Reinforcement Learning, Agent Technologies for Energy Systems[END_REF] and Huang et al. [START_REF] Huang | A Multiagent Minority-Game-Based Demand Response Management of Smart Buildings Toward Peak Load Reduction[END_REF] have designed MASs using Markov Decision Process (MDP) and game theory, respectively. A thorough study of existing literature suggests that most of the proposed MAS architectures can be classified as hierarchical, thus not fully decentralized. This property makes the system vulnerable to single point of failure (SPOF) issues and limits performance as sub-problems of limited size only, can be handled when strong real-time constraints are considered, which can be quite limiting for large-scale problems.

To overcome these challenges, a system based on a particular type of MAS, Adaptive Multi-Agent System (AMAS), is being developed. In AMAS, the goal is achieved through cooperative interactions of agents i.e., emergence. The main characteristics of such a system are autonomy, decentralization, sociability, and reactivity. The system is called ADaptive Energy Management In Smart Grids (ADEMIS) [START_REF] Blanc-Rouchosse | Multi-Agent System for Smart-Grid Control with Commitment Mismatch and Congestion[END_REF] and is capable of handling, by exploiting the storage capacity of EVs for ancillary services, basic grid stability challenges i.e., line congestions, as well as potential antagonist issues such as handling commitment mismatch while ensuring desired charging of the EVs.

It is equally important to keep a check of the obtained solutions' qualities, when developing a new methodology. For that purpose, the problem is simplified, modeled, and solved using Mixed Integer Linear Programming (MILP) optimization, which has already been used for optimal power flows in smart grids [START_REF] Liu | A MILP-based distribution optimal power flow model for microgrid operation[END_REF]. This enables to draw a comparison between a basic, uncoordinated charging strategy, the classical MILP optimization approach and the newly proposed AMAS methodology. A good comparison would allow us to better understand the position of AMAS solutions among different methodologies and quantify the room for improvements. Optimization formulation and linear, convex approximations for MILP along with AMAS modeling is presented in Section II. Following this, case studies are defined for testing in Section III and finally the results obtained are discussed in Section IV and conclusions are made.

II. SYSTEM MODELING

A. Optimization Formulation

The starting point is the power flow constraints of the network, which ensure power flow in the network based on its physical constraints. and represent the active and reactive powers at bus a at instant t. The power at each bus is equal to the sum of power injected or drawn from the bus.

, , and ∆ , stand for the day-ahead planned EV power connected to bus a at instant t, and the change in that power, hence a decision variable for optimization, to achieve the system's objective.

= - 1 = , + , - , 2 -! , , + ∆ , " = - 3 = , + , - , 4 -! , , + ∆ , "
Equation (5) links the power flows in the network to the bus voltages and network losses. % stands for the voltage at bus a, while & ' * represents the conjugate of the admittance matrix. Equations ( 6) and ( 7) relate powers at the bus with total inflows and outflows. Equations (8) and [START_REF] Sengupta | Oahu Solar Measurement Grid (1 YearArchive): 1-Second Solar Irradiance[END_REF] ensure that the power values remain within specified limits. Next comes the constraints on bus voltages and currents. These can be classified as Distribution System Operator (DSO) constraints. [START_REF]2030 climate & energy framework[END_REF] 11

% , ≤ |% | ≤ % , 1 10 |6 ' | ≤ 6 ',
EV constraints serve the purpose of modeling the dynamics of the EVs and also linking it with the network flow constraints. The State of Charge (SoC) of an EV is dependent on the previous state and charging or discharging power at the present interval. The SoC of the EV should be within specified limits, and it should be greater than the desired minimum SoC at the departure time. Similarly, the State of Health (SoH) of an EV, models health of the battery and should never be negative.

, 78 , and , 978 , is the EV charging and discharging power connected to bus a at instant t. Term : is a binary variable which makes sure that an EV is either charging or discharging, and not both, at interval t.

tan ∅ 20

The objective is to minimize the mismatch cost paid by the BRP, (21). BRP is charged for a negative mismatch (i.e., over-consumption) and vice versa. The negative mismatch cost is higher compared to the compensation for a positive mismatch. Terms is the real-time power consumption in the BRP perimeter, while is the dayahead planned power defined by (25). PV forecasts are usually provided as a one-hour average value for distributed generation, while for loads the resolution is smaller i.e., 10 minutes. The errors in the forecasts are also modeled here.

`ab ∆ >?,B = ^_ = 21 `ab ∆ >?,B ∑ *= 9 [ \ ] 9 [ \ -= 9 Q \ ] 9 Q \ + cdc eFJ cdL ] 9 [ \ = `fg !0, ∑ !* \ -^_ +∆ " d cQL ∆c[∆ dc∆c " 22 ] 9 Q \ = `fg !0, ∑ !* ^_ - \ +∆ " d cQL ∆c[∆ dc∆c " 23 ^_ = ∑ 24 \ = L ∆c ∑ `hfb L8 d cQL ∆c[∆ dc∆c * , 1 + hii + -`hfb LR , *1 + hii + -`hfb LR 1 + hii 25

B. MILP Approximations

To linearize the original non-linear and non-convex problem, ( 6) is approximated, obtained by applying small angle approximation ( sin l -l ' = l -l ' ), ( 26)-( 27). Also, the subtraction of the positive energy mismatch cost (a convex function), in the objective function makes the curvature of the total objective function non-convex. To retain the convexity, the positive mismatch cost has been considered as 0 €/MWh, which represents a worst-case.

' = m ' *% -% ' + + n ' *l -l ' + 26 ' = n ' *% -% ' + + m ' *l ' -l + 27 `ab ∆ >?,B = ^_ = `ab ∆ >?,B ∑ *= 9 [ \ ] 9 [ \ + cdc eFJ cdL

C. ADEMIS Modeling

This work is a continuation of Blanc-Rouchossé et al. [START_REF] Blanc-Rouchosse | Multi-Agent System for Smart-Grid Control with Commitment Mismatch and Congestion[END_REF]. The basic architecture of the system is the same as described in the previous work. Each physical component of the distribution grid is modeled as a cooperative agent, with each agent having its own objective, a partial information of the environment and a set of actions depending on the situation. The AMAS system executes in a loop. During each iteration, an agent calculates its own criticality, between 0 and 1, and compare it to the criticalities of the neighbouring agents. Subsequently it decides to take an action, either to satisfy its own objective or to help the neighbouring agents in achieving theirs. In the following paragraphs, the modelling of each type of agent is defined.

1) Line Agent: Line agents communicate criticalities, (29)-(30), based on the line current values. Term 6 8 is the threshold value at which line becomes critical and is set to 0.8 times the rated line current value 6 1 . =i , is the anticipated criticality of the line to ensure dynamics of the system, so that the line does not keep on oscillating between highly critical and non-critical states. Once a line becomes congested, its line agent will keep on communicating the anticipated criticality value so EV agents can continue cooperating with the line without disturbing the dynamics of the system.

3) EV Agent: EV agents calculate their respective criticalities based on (33). In (33),

1 is the maximum charging or discharging power of each EV,

, N O
is the time remaining before departure of the EV. After this calculation, they compare it to the criticality signals they have received from the neighbouring agents. Based on this, each EV agent makes the decision, either to help a line for congestion, or the BRP for commitment mismatch or to continue charging, by calculating the updated power for that interval. The updated power is calculated using (34). =i 1 stands for the maximum criticality received while =i 1,T is the second maximum criticality among antagonist (opposite) signals In the above equation, h is a function which tackles the situation when EV agents receives antagonistic requests. For example, if line is congested (thus requesting EVs to charge less) and BRP commitment mismatch is positive (thus creating favourable conditions for EVs to charge more) at a given instant simultaneously, then EV agents have to decide which request is more critical. Thus, the h function acts like a weighting function, given in (35).

ℎ =i , = ⎩ ⎪ ⎨ ⎪ ⎧ 1 + …[L c †,pEF =i , , as=i , ‡ \ , hgˆ! } B, † w B Qw G " , as \ , ≤ =i , ≤ \ , 1 1 - ‰ L[c †,pBq =i , , as=i , Š \ , 1 35
\ , 1 and \ , are thresholds for anticipated line criticality. Between two thresholds, the h parameter behaves exponentially and linearly otherwise. Terms t and t ' are two parameters dependent on values of ‹, OE, \ , 1 and \ , . The development of the methodology for ADEMIS being at an early stage, it can be observed that a simple strategy consisting in reducing the commitment mismatch was selected at first, thus not considering the potential reward for injecting power when there is an excess of PV production. However, the goal of this paper consists in assessing the room for improvement with such a usually envisaged strategy among BRPs, as well as to check that constraints are satisfied.

III. CASE STUDY

A. Distribution Network

The "IEEE Low Voltage Test Feeder" (IEEE LVTF) [8] is selected for this case study. As shown in Fig. 1, it consists of 3 zones (A, B &C). It is connected to a substation through a 11 kV/416 V transformer. The network consists of 55 households. These 55 household buses are also equipped with a PV and an EV connection, making the size of total EV fleet for this study equals to 55. For a congestion study, line 117 (highlighted in red) of the network is selected. For ADEMIS power flow, the network is simulated in DIgSILENT PowerFactory (PF).

B. Datasets

This study is carried out for a total time of one hour (4 imbalance settlement periods of 15 minutes each), from 16:00 till 17:00. The dataset for load profiles of 55 households is provided along with the IEEE LVTF model. The resolution of the data is 1 minute. The irradiance data for PV profiles (resolution of 1 second) is obtained from the database of the National Renewable Energy Laboratory (NREL) [START_REF] Sengupta | Oahu Solar Measurement Grid (1 YearArchive): 1-Second Solar Irradiance[END_REF]. The power generated by a panel is calculated using (36). 6ii is the irradiance at time instant t, A is the area of the panel, 10 `T here, I , the efficiency of the panel, is equal to 0.17. The error for day-ahead PV forecast is set to 15% during the studied hour. The EV fleet modeled for this system can be classified into three distinct types of EVs, 1) EVs arriving at start of period 1 for charging and departing soon, which would create a congestion in the network line 117, 2) EVs arriving during the simulation time and departing near the end or after simulation time, and 3) EVs connected to the network before the start of the simulation and departing afterwards. All EVs present the same capacity of 30 kWh and have a maximum charging or discharging power of 7 kW. Terms ;Ž= and ;Ž= 1 have been set to 0.3 and 0.8 respectively, while ;Ž= * N O + is equal to 0.7.

= •I 6ii 36

Additional data for imbalance costs of the BRP is obtained through the French TSO Réseau de Transport d'Électricité (RTE) [START_REF]Réseau de Transport d'Électricité (RTE) balancing costs[END_REF]. The data provided by the TSO has a resolution of 30 minutes, thus a linear resampling of data is done to obtain costs values for each 15 minutes period.

C. Antagonist Issues

This study is modeled to incorporate a variety of scenarios to make it more realistic and non-trivial. During period 1, PVs are over-producing, while under-producing in the remainder. BRP mismatch is present in all 4 periods, due to differences in the forecasted and real-time values. Also, during period 1 a line congestion is expected when uncoordinated charging occurs. The rated value for line current is 0.065 kA. The line type (number of cores and cross-sectional area) is also given with the IEEE LVTF model. The rated current value is obtained from the technical sheet [START_REF]Low voltage power cables[END_REF] for the congested line. To provoke a congestion in the network, 8 EVs in zone A and 4 EVs in zone B are simulated to arrive near the start of the simulation. A summary is presented in Fig. 2. EVs causing congestion are charging at their maximum specified power as their departure time is around 16:12. As the EVs load forecast is made by taking the average over 10 minutes-periods, if the total load is very high (enough to cause congestion) for a portion of 10 minutes and very low for the remainder, the average value would be in the middle and would not depict congestion for the portion when a very high total load is present, which is exactly what happening in this scenario. As stated, we have over-production and congestion in period 1, which gives rise to an antagonistic scenario, where EVs may be incentivized to consume more power but as there is a line congestion, line agent would require EVs to draw lesser power through the line.

IV. RESULTS

In this section, first the results obtained from MILP and ADEMIS are discussed individually and then a comparison of both along with the basic charging strategy is made.

A. MILP Results

To model MILP formulation CVXPY [START_REF] Diamond | CVXPY: A Python-embedded modeling language for convex optimization[END_REF], with GUROBI solver, is used. The computing machine is equipped with Intel i5-10210U CPU with base speed of 2.11 GHz and 8 GB RAM (DDR4). The temporal resolution of the study is 1 minute. Fig. 4-b shows that during period 1, when congestion is expected to occur with an uncoordinated charging, the total EV consumption is lower than the forecasted consumption. Power flow from the grid to the congested zones is reduced, and the EVs present in both zones, which do not need to depart early, started supporting the EVs which are departing soon by charging less or discharging.

To check the impact of the approximations on the quality of the solution, the optimized power profiles obtained through MILP are simulated in the IEEE LVTF PowerFactory model. Active power grid profiles are shown in Fig. 3-a. The percentage cosine similarity between both curves is 99.018%. This confirms a minor impact of MILP approximations on the quality of the solution. To draw a fair comparison, battery aging cost has been omitted in the MILP formulation, as it is not present in the current version of ADEMIS. This results in a significant number of oscillations in the solution profiles.

B. ADEMIS Results

Agents' behaviour is implemented in JAVA. Power flow for each iteration is done in PowerFactory. Both platforms communicate through a python script. As the agents react to the system in real-time, the time step considered for each iteration is 1 second. The system consists of 55 EV agents. Fig. 4-a confirms that the congestion in the network is also avoided for ADEMIS solution as the line current is smaller than the rated current value of the line. Present version of ADEMIS is configured to reduce the commitment mismatch, as explained in (31-32). Further, the EVs present in both zones, which do not need to depart soon, are helping the EVs departing soon by discharging.

C. Results Comparison

Fig. 4-a shows the line currents in line 117. Both ADEMIS and MILP solutions are successful in managing the line congestion, which should have occurred with a basic charging strategy. Further, both for MILP and ADEMIS, the SoC criteria (minimum SoC reached before EV departure) of all the EVs are satisfied, which cannot be guaranteed for the basic charging approach in case a zone(s) gets disconnected from the grid side due to congestion or other failures in the network.

A comparison of energy mismatch for each period is presented in Fig. 5. This comparison is between different formulations presented earlier, obtained through different methodologies to tackle the same set of problems. Negative values indicate positive energy mismatch (over-production), whereas positive values mean negative energy mismatch (over-consumption). For periods 2,3 and 4, ADEMIS and MILP follow the same trend of significantly reducing the expected negative energy mismatch, compared to the uncoordinated charging strategy. In period 1, the cost of injecting more power and battery aging is zero in this preliminary study. Hence, the MILP approach proposes one feasible solution among others, where a large amount of power is injected in the grid (at no cost). However, using the non-convex objective function, this injection is rewarded. Under these conditions, the behaviour of ADEMIS is counterproductive compared to the basic charging, as it reduces the injection of power into the external grid. The goal of studying the conditions in period 1 was primarily to check that the EVs and DSO constraints were satisfied when antagonistic requests were received by the agents, knowing that the basic strategy implemented in the agents would not be relevant in this particular case. More sophisticated BRP functions are considered for future work in order to improve the performance of the overall system.

Finally, the cost comparison for each formulation is done by placing the obtained solution values in (21) (for MILP both (21) and (28) are used) and then calculating ratios with MILP PF results as reference. Again, it can be verified, in Table I., that the solution obtained through ADEMIS stands between MILP and basic charging strategy. MILP, (28), solution is closer to the ADEMIS solution, as in neither of the approaches the exploitation of the reward by injecting power is sought.

V. CONCLUSION

A comparison of MILP optimization and AMAS has been presented to solve the grid stability and commitment mismatch problems in the distribution grid. Both approaches efficiently handle the given challenges. Compared to the basic charging strategy, ADEMIS solution is performing better overall, but a substantial gap between ADEMIS and MILP solutions indicates that there is large room for improvement. On the flip side, AMAS is more modular and scalable compared to the MILP. Anticipative stochastic optimization is based on forecasts which can have errors and requires a large computing time, being NP-complete. Unlike this approach, an AMAS system reacts to the environment in real-time and the computing time of each agent is independent of total number of agents in the system. This provides motivation for the future work, which includes improving the behaviour of the agents as well as bus voltage management to study large-scale unbalanced networks. This would give a better benchmark to compare not only the quality of the solutions but also the time complexity to obtain those solutions.

  and positive mismatch costs and energy values, respectively. ^_

  Agent: BRP agent sends its criticality, (31)-(32), based on its instantaneous mismatch value. Terms t and t ^_ are tuning parameters with arbitrarily selected values 0.01 and 100, respectively.

Fig. 1 .

 1 Fig. 1. Graphical representation of the IEEE low voltage test feeder

Fig. 2 .

 2 Fig. 2. Summary of the case study, EVs arriving (↑) and departing (↓)

Fig. 3 .

 3 Fig. 3. MILP and PowerFactory results comparison

Fig. 4 .

 4 Fig. 4. Results comparison: a) Line 117 current b) Total EV powers

Fig. 5 .

 5 Fig. 5. Total energy mismatch comparison for each period

  Negative mismatch cost values are 91.26, 108.83, 103.28 and 92.124 €/MWh, while positive mismatch values are 82.56, 98.46, 93.44, 89.35 €/MWh for periods 1, 2, 3 and 4, respectively.

TABLE I .

 I COST COMPARISON (RATIO) OF STUDIED FORMULATIONS

	Ratio			Strategy
	• • ''" •	"	Basic	ADEMIS	MILP (28)	MILP (21)	MILP (PF)
	Period 1	0.165	0.063	0	1.006	1
	Period 2	-0.933	0.733	0	0	1
	Period 3	-37	-22	0	0	1
	Period 4	-36	-8	0	0	1
	Total		-0.21	-0.02	0	0.95	1

ACKNOWLEDGMENT

This work is carried out as part of the "Opt-Real" project funded by the Brittany Regional Council which is gratefully acknowledged.


N Massacret

J Moysan

M A Ploix

J P Jeannot

G Corneloup

Modelling of ultrasonic propagation in turbulent liquid sodium with temperature gradient

The use of ultrasonic instrumentation in sodium-cooled fast reactors (SFR) requires to understand and to predict how ultrasonic waves can be deflected, slowed down or speeded up, depending on the thermo-hydraulic characteristics of the liquid sodium. These thermo-hydraulic characteristics are mainly the local temperature and flow speed of the sodium. In this study we show that ray theory can be used to simulate ultrasonic propagation in a medium similar to the core of a sodium-cooled fast reactor, in order to study ultrasonic instrumentation and prepare it installation and utilisation in the sodium of the nuclear reactor. A suitable model has been developed and a set of thermo-hydraulics data has been created, taking account of the particularities of the sodium flow. The results of these simulations are then analysed within the framework of acoustic thermometry, in order to determine which disturbance must be taken into account for the correct operation of the temperature measurement.

 . It is possible to obtain excellent measurement dynamics: as the speed of ultrasonic waves is around 2300

I. INTRODUCTION

Sodium-cooled fast reactors have been chosen by France as the reference option in the framework of the Generation IV international forum. These reactors must achieve the Generation IV objectives in terms of sustainable development, economic competitiveness, reliability, proliferation resistance and protection of people 1 .

The CEA (atomic energy and alternative energies commission), EDF (Electricité de France) and AREVA have formed a partnership to coordinate their research programmes and to develop a fourth generation nuclear reactor prototype.

One of the many developments necessary to reach the industrial standard is a research programme named ISIR (In-Service Inspection and Repair) which has been initiated by the three partners 2 . One of the objectives of this research programme is to take the configuration of this new reactor and the characteristics of the coolant into account in order to develop suitable instrumentation for continuous in-service monitoring. Actually in these reactors, the coolant is liquid sodium which is an opaque liquid metal, preventing the use of conventional optical methods and reaching temperatures around 550°C in normal operating. Various instrumentation techniques a : Author to whom correspondence should be addressed. Electronic mail: joseph.moysan@univ-amu.fr associated with the use of sodium are studied in this programme, including detection of bubbles in the sodium, ultrasound transmission through interfaces, and acoustic thermometry [START_REF] Cavaro | Proceedings of Meetings on Acoustics[END_REF][START_REF] Paumel | [END_REF][START_REF] Massacret | AIP Conference Proceedings[END_REF] .

Acoustic thermometry is a technique based on the temperature-dependent speed of ultrasonic waves in a medium. It is possible to propagate the ultrasounds in a solid wave-guide immersed in the hot fluid, as in the case of some acoustic thermometers [START_REF] Daw | Proc. NPICHMIT 2010[END_REF] , or propagate them directly in the fluid: it is this second application that is under consideration in this paper. High temperature transducers are also studied [START_REF] Baba | [END_REF]8,[START_REF] Parks | Rev. of Progress in QNDE, AIP Conf. Proc[END_REF] . We plan to develop a method involving the propagation of an ultrasonic wave to two surfaces, separated by known distance, which generate echoes. The mean local speed of the medium is obtained by measuring the ultrasonic wave time of flying between these two surfaces. As the relationship between the speed and the temperature is known, the temperature of the sodium is obtained. Leibowitz et al. have established a relationship between ultrasound speed and temperature in sodium [START_REF] Leibowit | [END_REF] :

𝑐 = 2577.3 -0.536 ⋅ 𝑇 ( 1
)
where 𝑐 is the speed of the ultrasonic waves in meters per second and 𝑇 is the temperature of the sodium in degrees Celsius. This formula is valid between 182°C and 896°C with an accuracy of 0.4% in pure sodium.

For in-reactor application, the edges of the fuel subassembly heads are good candidates for generating the echoes. It would then be possible to measure the temperature of the sodium directly on the diameter of the subassembly outlets. An example of an experimental configuration is shown in figure 1.

FIG. 1. Acoustic thermometry principle with a grazing incidence beam.

m.s -1 in sodium at 550°C, the response time is mainly due to the signal processing time. Moreover it would be possible, with only one transducer, to measure the temperature of the sodium exiting several subassemblies. This was proved in water experiment, where the temperature of the fluid exiting a row of six subassemblies was measured, using a single transducer emitting a grazing incidence beam [START_REF] Taylor | AEA Technology[END_REF] . This method would significantly reduce the amount of instrumentation above the core.

However the thermo-hydraulic conditions above the core may disturb the propagation of the ultrasonic waves between the ultrasonic transducer and the subassembly heads, in terms of time delay and deflection.

There are several sources of thermal inhomogeneities above the core: the temperature difference of the sodium flowing out two neighbouring subassemblies could reach 50°C following the core conception. Moreover the sodium circulating in the spaces between the subassemblies, and the one flowing out the spaces left clear for insertion of control rods or safety devices, are cooler by several tens of degrees Celsius, than the sodium flowing out the subassemblies. The ultrasonic waves will therefore propagate in a medium in which the temperature is very inhomogeneous.

In addition, the flow above the core is turbulent with local speeds reaching approximately 3 m.s -1 , and speed gradients are about several metres per second per centimetre. This turbulent field may also have consequences on the propagation of the ultrasonic waves. This phenomena is used in acoustic debitmeter to measure the speed flow [START_REF] Liu | [END_REF][14] . In the case of acoustic thermometry, it could lead to errors in the measurement of the temperature.

In order to be able to anticipate and deal with these potential problems, a 2D ultrasonic propagation model, suitable for this medium and this instrumentation, has been studied to simulate the path of the ultrasonic waves.

These simulations will enable the disturbance caused by the thermo-hydraulic characteristics of sodium to be quantified, and will determine of whether they are problematic for acoustic thermometry. In this article, the reasons for choosing this model will be explained, then its theory detailed. Finally, various simulations highlighting probable types of disturbance and their effects on acoustic thermometry will be analysed.

II. THE ACOUSTIC MODEL

A. Choice of the appropriate model

It is important to take the thermo-hydraulic characteristics of the medium into account when choosing an appropriate, effective ultrasonic propagation model. Some codes have been developed for extremely turbulent media in which the Mach numbers are close to one, using methods involving different degrees of parabolic approximation 15 . Despite the precision of these codes, the sodium flow above the core of an SFR is not fast enough in comparison with the speed of ultrasonic waves to justify the use of such complex methods. For this study, in which the Mach number is around 10 -3 , it is possible to use the frozen fluids hypothesis, in which the thermohydraulic conditions of the medium are considered to be stationary when an ultrasonic wave is propagated. This approximation is used in particular for ultrasonic studies of turbulence in fluid media 16 .

Another important parameter is the ratio between the characteristic size of the inhomogeneities of the medium and the ultrasonic wavelength. The size of the thermal inhomogeneities depends mainly on the thermal conductivity of sodium. As it is a liquid metal, its thermal conductivity is high and its Prandtl number is very low, around 10 -3 . The Prandtl number is the ratio between kinematic viscosity and thermal diffusivity. It is used to determine whether the temperature field of a fluid will be significantly influenced by its speed field (high Prandtl number), or influenced by it very little (low Prandtl number). Small thermal inhomogeneities will therefore disappear quickly in this medium. Several experiments have been carried out in sodium, and it seems that the characteristic size of the thermal inhomogeneities is more than 2 centimetres [17][START_REF] Moro | IAEA International Working Group on Fast Reactors, Specialist Meeting on Correlation Between Mateial Properties and Thermo-Hydraulics Conditions in LMFBRs[END_REF] . The velocity of a 2.25 MHz ultrasonic wave propagating in sodium at 550°C being around 2300 m.s -1 , its wavelength is around a millimetre. It is therefore possible to make the assumption of high frequency waves 14 :

𝜆 ≪ 𝐿 and √𝜆𝑋 ≪ 𝐿 (2)
where 𝜆 is the ultrasonic wavelength, 𝐿 the characteristic size of the thermal inhomogeneities, and 𝑋 the length of the path travelled by the ultrasonic waves in the fluid.

With these two assumptions (frozen fluid and high frequency), the acoustic ray model based on geometric acoustic theory must be favoured 14 .

Another important parameter which encourages the choice of this model is the calculation cost. A large number of simulation codes use the finite or spectral element method [START_REF] Zhang | [END_REF][20][21][START_REF] Moussou | Flow Induced Vibration -International Conference[END_REF] . These codes are very efficient for simulating the interference of an ultrasonic wave with inhomogeneities and the distribution of the acoustic amplitude in the medium. However, calculating the propagation of a wave of several megahertz over several tens of centimetres generates significant calculation costs.

The application of the ray tracing method request generally few numerical resources because the calculations concern only the small part of the geometry in the area of the propagating ray. Moreover the calculation of the path of each ray is independent, so it is possible to calculate these paths in parallel on several processors. In order to determine the distribution of the acoustic amplitude, it is then possible to apply the Gaussian beam summation method.

B. Acoustic ray model

The geometric acoustic ray model consists on following the spatial progress of a point on a wavefront. The path travelled by this point is called the acoustic ray. A wavefront is defined by all the particles in the propagation medium which, at a given time, are set in motion by the passing of an alternation. An iterative method is used in this study: at each step of the calculation, the local thermo-hydraulic characteristics (temperature and flow speed of the sodium) are used to calculate the new direction of propagation of the wave. The path of the ray is then calculated according to this direction over a given distance. The same operation is carried out at this new position.

Thus, the two most important steps at each iteration, are taking the thermo-hydraulic data into account and calculating the direction of propagation.

The calculation model is based on the fact that if a particle in the medium is moving at a flow speed 𝑣 ⃗(𝑥 ⃗), the wave will propagate at a speed 𝑣 ⃗ + 𝑐 ⋅ 𝑁 ⃗ , where 𝑁 ⃗ is the local direction of propagation of the wave (i.e. normal to the wavefront) and 𝑐 is the velocity of the ultrasonic waves in an equivalent static fluid [START_REF] Pierce | Acoustical Soc of America[END_REF] . With 𝑥 ⃗ (𝑡) a point on the wavefront, this point moves at speed 𝑣 ⃗ + 𝑐. 𝑁 ⃗ , thus:

𝑑𝑥 ⃗ 𝑑𝑡 = 𝑣 ⃗ 𝑥 ⃗ , 𝑡 + 𝑁 ⃗ 𝑥 ⃗ , 𝑡 . 𝑐(𝑥 ⃗ , 𝑡) (3)
Then let 𝜏 is the function describing the coordinates of each point of the wavefront. At instant 𝑡, all the points 𝑥 ⃗ belong to the wavefront and are described by this function. Thus:

𝑡 = 𝜏 𝑥 ⃗ (4)
If the time interval Δ𝑡 between two positions on the wavefront is short enough, it is possible to write:

𝑡 + Δ𝑡 ≈ 𝜏 𝑥 ⃗ (𝑡) + 𝑥 ⃗ ̇ (𝑡)Δ𝑡 ≈ 𝜏 𝑥 ⃗ + Δt 𝑥 ⃗ ̇ (𝑡) ⋅ ∇ ⃗ 𝜏 (5)
I.e. ∇ ⃗ 𝜏 = 𝑠 ⃗, with 𝑠 ⃗ is the slowness vector:

𝑠 ⃗ = 𝑁 ⃗ 𝑐 + 𝑣 ⃗ ⋅ 𝑁 ⃗ (6)
Since 𝑠 ⃗ is dependent on 𝑥 ⃗ (𝑡):

𝑑𝑠 ⃗ 𝑥 ⃗ 𝑑𝑡 = ∇ ⃗ 𝑠 ⃗ ⋅ 𝑥 ⃗ ̇ = 𝑐 𝑁 ⃗ ⋅ ∇ ⃗ 𝑠 ⃗ + 𝑣 ⃗ ⋅ ∇ ⃗ 𝑠 ⃗ (7)
And:

𝑑𝑠 ⃗ 𝑥 ⃗ 𝑑𝑡 = 𝑐 Ω 𝑠 ⃗ ⋅ ∇ ⃗ 𝑠 ⃗ + 𝑣 ⃗ ⋅ ∇ ⃗ 𝑠 ⃗ With Ω = ⃗⋅ ⃗ = 1 -𝑣 ⃗. 𝑠 ⃗ (8)
It is possible to define the location of any point 𝑥 ⃗ (𝑡) on the wavefront with a system using two equations:

⎩ ⎪ ⎨ ⎪ ⎧ 𝑑𝑥 ⃗ 𝑑𝑡 = 𝑣 ⃗ 𝑥 ⃗ , 𝑡 + 𝑐 𝑥 ⃗ , 𝑡 Ω 𝑠 ⃗ 𝑑𝑠 ⃗ 𝑥 ⃗ 𝑑𝑡 = - Ω 𝑐 ∇ ⃗ 𝑐 -𝑠 ⃗ ∧ ∇ ⃗ ∧ 𝑣 ⃗ -𝑠 ⃗ ⋅ ∇ ⃗ 𝑣 ⃗ (9)
Thus, using these two equations (9), the path of an acoustic ray can be determined.

C. Gaussian beam summation

Ray theory cannot be used on its own to determine the amplitude of the acoustic signal at all the points of the geometry. This amplitude is considered to be zero out of the rays, which is not the physical reality.

To overcome this drawback the Gaussian beam summation method can be applied 14,[START_REF] Cerveny | [END_REF][25][26] . The first step in this method involves calculating the propagation of a beam of acoustic rays. Each beam is then given a Gaussian amplitude profile. This profile is characterised by a Gaussian decay according to the distance from the ray. Its width at mid-height and its maximum amplitude are two parameters which change along the ray according to the frequency of the propagated wave. Then, at the point at which it is necessary to know the acoustic amplitude of the signal, the amplitude from each ray is summed. With this Gaussian decay of the energy, the further a ray is from the point of interest, the weaker its energy contribution at this point.

Again, the advantage of this method is its speed: it is not necessary to apply it along the whole length of the acoustic beam to obtain the amplitude at a specific point; it simply needs to be applied in the area concerned and at a chosen moment.

The starting point of the calculation is the conventional wave equation, here in a 2-dimensional medium [START_REF] Cerveny | [END_REF] :

∂ u ∂x + ∂ u ∂y = 1 c (∂ u) ∂t ( 10
)
where 𝑢 is a physical parameter linked to the movement of the wave, 𝑐(𝑥,y) is the speed of the wave and 𝑡 is the time.

The first step in the demonstration is to change the system of coordinates: the Cartesian coordinates system (here 𝑥 and 𝑦) is abandoned in favour of the coordinates linked to the ray in question. Giving 𝑡 ⃗ the direction vector which is tangent to the ray, and 𝑛 ⃗ the direction vector which is orthogonal to it. The scheme in FIG 2 . Coordinate system linked to the ray.

explains this change of coordinates.

. The scalar 𝑠 also appears in this diagram. This indicates the curvilinear abscissa of the ray. 𝑠 is incremented at each iteration when the ray method is applied. Point A is a point of interest in the medium and A' is the position of point A following an iteration. In the coordinates system (𝑠,n) the wave equation ( 10) becomes:

1 ℎ 𝑢 , + ℎ 𝑢 , - ℎ 𝑐 𝑢 , + 𝑢 , 1 ℎ + 𝑢 , ℎ , = 0 (11)
With:

ℎ = ( , ) ( , )
For formula (11), the partial derivatives are indicated using a comma. This notation will be used in the remainder of this article.

Then the parabolic approximation method is applied. This method enables only those waves propagating in a small aperture cone centred on the axis of the ray to be taken into account in the calculation. Waves created during reflection or diffraction are thus ignored. This method reduces the calculation cost substantially. It is also fully justifiable here, as this study is only examining the propagation of the waves: interactions with the solid structures in the medium are not being studied. In addition, since the assumption of high frequency waves is applied, the medium in which the calculations are carried out does not contain any thermo-hydraulic phenomena which could cause diffraction or diffusion.

The application of this method consists of stating 𝑢 in the form:

𝑢(𝑠, 𝑛, 𝑡) = 𝑈(𝑠, 𝑛, 𝑤). 𝑒𝑥𝑝 -𝑖𝑤 𝑡 - 𝑑𝑠 𝑐(𝑠) (12)
Where 𝑡 is the time, 𝑤 is the angular frequency of the signal and 𝑈 is a constant of the same type as 𝑢, specific to each ray.

A condition is thus placed on 𝑛 in the steps of the calculation:

𝑛 = 0 𝑤 (13)
Where the value ½ has been defined empirically [START_REF] Babič | Alpha Science International[END_REF] . The notation 0 means that 𝑛 wouldn't be much more superior to 𝑤 and therefore that the scope of the method is only suitable for an area close to the ray, which decreases as the angular frequency increases.

Following the resolution steps of Cerveny 24 , a system of two linear differential equations is introduced by inserting two complex functions 𝑝(𝑠) and 𝑞(𝑠), and a real value Φ which can differ from one ray to another. The amplitude of the acoustic signal 𝑢(𝑠,n,t) at a certain distance 𝑛 from the ray is thus calculated using a three-equation system:

⎩ ⎪ ⎨ ⎪ ⎧ 𝑞 , = 𝑐𝑝 𝑝 , = - 𝑞 𝑐 𝑐 , 𝑢(𝑠, 𝑛, 𝑡) = 𝛷 𝑐(𝑠) 𝑞(𝑠) 𝑒𝑥𝑝 -𝑖𝜔 𝑡 - 𝑑𝑠 𝑐(𝑠) + 𝑖𝜔 2
inhomogeneities (temperature or speed) which occur in this medium are around a centimetre and there would therefore be significant redundancy of data 29 .

It was therefore decided to carry out a linear interpolation of the thermo-hydraulic data to determine the temperature and speed conditions in the location through which the ray passes. In this way, it is possible to use very varied sets of thermo-hydraulic data (with fine or coarse grids depending on the turbulence of the medium), while maintaining temperature or flow speed variations descriptions that are fine enough to allow propagation of the ray.

Moreover it is possible to use results of thermo-hydraulics simulations with various meshing (structured or not).

This is an interesting point for research gathering thermo-hydraulics and acoustics issues.

In the following first simulations given in this article, digitally generated sets of thermo-hydraulic data will be used. Thus we propose a simplified model which reproduced the reality in terms of size and amplitude of the inhomogeneities.

B. Simulated media

As shown previously, the thermo-hydraulic conditions of the medium will influence the propagation of ultrasonic waves in liquid sodium. It is therefore important to know to what extent this disturbance will occur and if it will affect ultrasonic propagation significantly. Three simple cases are simulated in this article, with the aim of studying the temperature and flow effects separately. The medium is the liquid sodium just above the subassembly heads.

When carrying out acoustic measurement of the temperature of the sodium at the outlet of a subassembly, an ultrasonic beam propagates above the core until it is reflected on the edges of the subassembly. This beam successively encounters hot areas and cold areas, created respectively by the sodium exiting the subassemblies and the spaces between the subassemblies. The shape of the sodium thermal and flow inhomogeneities will disturb the ultrasonic waves to a greater or lesser extent depending on the direction of ultrasonic propagation. This medium also has an inhomogeneous flow speed field map: the hot sodium exiting the subassemblies moves at a speed of around 3 m.s -1 , creating a significant speed gradient with the sodium moving between the subassemblies at a much lower speed (close to 0 m.s -1 ).

The first medium simulated in this study is an horizontal cross-section of a row of three subassemblies, discharging hot sodium at 570 °C. The temperature of the sodium in the spaces between the subassemblies is 500°C. The thermal gradient is linear and extends over a distance equal to the thickness of one subassembly edge.

A map of the temperatures of this geometry, is shown in figure 4. The flow speed is not taken into account in this case, only the effect of the temperature. The aim here is to highlight the influence of the specific nature of the temperature field just above the subassemblies: in this situation there is no random distribution of the temperature which changes over time, but relatively static temperature gradients. It is thus interesting to observe the accumulation of these organised and repeated inhomogeneities on an ultrasonic beam.

The second geometry used for the simulations will additionally take into account certain subassemblies adjacent to a row of subassemblies on which ultrasonic thermometry is being carried out. The geometry of the second case is shown in FIG 5. As previously, the sodium exiting the subassembly heads is at 570°C and that in the spaces between the subassemblies is at 500°C. The thermal inhomogeneities are stable over time and space and the ultrasonic beam will therefore be propagated through a grid of thermal inhomogeneities. Here also the movement of the sodium is not taken into account. The aim of the corresponding simulation is to determine whether a beam propagating through this grid will adopt specific directions of propagation, which could disturb the thermometry of certain subassembly heads.

The third case studied is a vertical cross-section of the speed field of the sodium from three aligned subassemblies. This map is shown in FIG 6. This geometry highlights the three ascending vertical flows of sodium exiting each of the three subassemblies.

The speed of this flow is 3 m.s -1 . The sodium between the subassemblies is considered to be stationary. In the case of acoustic thermometry, an ultrasonic beam propagates through these flows up to the subassembly heads under consideration. The simulation to be carried out with this geometry does not take into account the effect of the temperature variations nor the turbulent movements which may occur in the mixed layer between the moving and stationary fluids. Its aim is simply to observe the cumulative effect of the repeated strong speed gradients on the path of beams with several incidence angles, and to see whether the effect is problematic for thermometry. The conclusions of these simulations will enable us to determine whether it is necessary to add this effect to that of the temperature in the cases covered in this study, or whether the effect of the temperature alone is adequate to describe wave propagation.

IV. RESULTS AND DISCUSSIONS

The first simulation below consists in ultrasonic rays propagating along a row of subassemblies (case 1). One hundred parallel rays are generated at an abscissa close to zero to obtain a beam in the axis of the row of subassemblies. The diameter of the beam is 50 mm. The results of this simulation are given in FIG 7.A) and B). .A shows the path followed by the simulated rays. Ultrasonic deflections appear in the areas in which propagating rays are not collinear with the gradient. This simulation shows that the rays curve towards the centre of the row of subassemblies. This is due to the fact that the ultrasonic waves are deflected towards the areas in which the speed of the medium is the lowest 26 . Here the sodium is hotter at the centre of the subassemblies and according to equation (1), it is in this area that ultrasonic waves propagate most slowly. heads, it is unlikely that echoes returning directly towards the transducer will be generated. The consequence of this orientation error could therefore be the loss of the ultrasonic signals used for measuring the temperature. For the same reason it is unlikely that echoes from rays propagating in the second row will be obtained. This simulation also indicates the importance of propagating an ultrasonic beam in the central axis above a row of subassemblies to reach a subassembly head and carry out a temperature measurement. Targeting for a subassembly head without taking this condition into account with another acoustic path, or using a transducer emitting an ultrasonic beam with too large aperture, would significantly complicate the processing of the received signal.

With the third set of thermo-hydraulic data, the effect of the flow speed on the propagation of the ultrasonic waves is studied. In these simulations five rays are propagated at angles of 5 ° , 2.5 ° , 0 ° , -2.5 ° and -5 ° in relation to the axis 0x. These rays appeared in figure FIG 9A.

A first simulation is performed in which the speed of the flow is not taken into account, and the five 5 rays meet at a point. These rays appear as crosses in figure FIG 9B. This is the centre of the arc of circle from which the rays are emitted. In this way all the rays have travelled the same distance. Furthermore the time of fly of each ray to reach the focus point is also different. The rays propagating with a downward angle relative to the axis Ox was slowed, and rays propagating with a bottom angle were accelerated.

These results appear in the These changes in the time of flight of rays come from the component of the flow velocity which is added or subtracted to the speed of ultrasound.

In the case of acoustic thermometry, where the ultrasonic beam has a downward trajectory to the sub-assembly head, then is reflected and returns with an upward trajectory to the transducer, changes in time of flight will be subtracted [START_REF] Leibowit | [END_REF] .

Moreover considering that the vertical deflection undergone by these rays is less than one millimetre, it seems that the influence of the flow velocity, in these specific thermo-hydraulic and acoustic conditions, are too minor to have a significant impact on the ultrasonic thermometry. It will however be necessary to take account of other turbulent hydraulic phenomena that exist above the core (impact of the core cover plug, turbulence near the control and command rods, etc.), before generalising this conclusion.

V. CONCLUSION

FIG 2 .

 2 FIG 2. Coordinate system linked to the ray.

FIG 4 .

 4 FIG 4. Thermal map of three subassemblies (Case 1).

FIG 5 .

 5 FIG 5. Thermal map of a grid of subassemblies (Case 2).

FIG 6 .

 6 FIG 6. Ascendant speed map of three subassemblies -Vertical cross-section (Case 3).

FIG 7 .

 7 FIG 7. A) Propagation of rays along a row of subassemblies. B) Wavefront obtained with the Gaussian beam summation method. FIG 7.A shows the path followed by the simulated rays. Ultrasonic deflections appear in the areas in which

FIG 7 .

 7 FIG 7.B shows the amplitude of the wavefront at a given time calculated with the Gaussian beam summation method. This figure presents the focusing of the beam in axis of the row of subassemblies and makes it possible to see whether the propagation of the rays in various areas in which the speed varies, causes a delay on certain parts of the wavefront. The wavefront shown inFIG 7.B has a slightly concave shape due to the faster propagation of the rays on each side of the beam, but this deformation remains weak. Thus, depending on the shape of the echoes

FIG 9 .

 9 FIG 9.A) Ray propagation across ascendant flows. B) Zoom on the focusing point of the rays. Crosses: calculation of the ray path in a static medium. Line: calculation of the ray path takes into account the influence of the flow velocity.

TABLE 1 .TABLE 1 .

 11 Variation in time of flight to reach the focusing point.

	𝑃𝑟𝑜𝑝𝑎𝑔𝑎𝑡𝑖𝑜𝑛 𝑎𝑛𝑔𝑙𝑒	𝑉𝑎𝑟𝑖𝑎𝑡𝑖𝑜𝑛 𝑖𝑛 𝑡𝑖𝑚𝑒 𝑜𝑓 𝑓𝑙𝑖𝑔ℎ𝑡
	𝑖𝑛 𝑟𝑒𝑙𝑎𝑡𝑖𝑜𝑛 𝑡𝑜 𝑎𝑥𝑒 𝑂𝑥 (𝑑𝑒𝑔𝑟𝑒𝑒𝑠)	𝑡𝑜 𝑟𝑒𝑎𝑐ℎ 𝑡ℎ𝑒 𝑓𝑜𝑐𝑢𝑠𝑖𝑛𝑔 𝑝𝑜𝑖𝑛𝑡 (𝑠)
	-5	1.40 ⋅ 10
	-2.5	6.93 ⋅ 10
	0	-9.0 ⋅ 10
	2.5	-7.13 ⋅ 10
	5	-1.41 ⋅ 10

Author Pre-Print -Journal of Applied Physics -May 2014Journal of Applied Physics 115(20):204905-204905-8 DOI: 10.1063/1.4875876

Finally the contributions of each ray at a particular point in the simulated medium must be summed in order to ascertain the acoustic energy at this point. In order to be able to observe the changes in a wavefront, or carry out a simulation with an acoustic signal as close as possible to reality, this code has been implemented so as to give Φ a real value that can be time-dependent. For the following simulations, in order to simplify the results, a very simple signal is used: several alternations of 2.

III. THERMO-HYDRAULIC DATA FOR THE SIMULATED MEDIUM A. Nature of the data

To preserve the speed of the code, it has been necessary to handle the thermo-hydraulic data in a specific way.

The frequency of the waves propagated in this code is in the region of a few megahertz. To simulate their propagation correctly, it is recommended that for each calculating step a ray is propagated on a distance of , i.e.

a distance close to a tenth of a millimetre [START_REF] Karweit | [END_REF] . In addition, as the path of a ray is sensitive to variations in the thermohydraulic conditions of the medium, it is not possible to propagate it in a coarse set of thermo-hydraulic data in which these conditions vary suddenly from one area to another with the risk of no longer complying with the assumption of high frequency waves (see Section 2.A.).

It is possible to use a set of thermo-hydraulic data containing grids close to a tenth of a millimetre, but this will generate a huge amount of data and lots of these will be unnecessary. In fact the characteristic lengths of the In this study a simulation code has been implemented in order to calculate the path follow by an acoustic wave in a turbulent sodium flow and inhomogeneous temperature. The application of the method of Gaussian beams summation also enables to determine the amplitude of the wavefront of the acoustic beam at various time of propagation.

The establishment of SFR core instrumentation such as acoustic thermometry is long and complex. The results calculated by this code are thus very useful to quickly test several positions or orientations of the transducer, or assess the impact of various thermal-hydraulic parameters on the temperature measurement. The first results presented in this paper indicate that it is preferable to orientate the transducer in the axis of a row of assemblies in order to achieve a particular assembly. These results predict concentration of ultrasounds energy along the axis of that row, and thus in the area where the echoes on the edges of assemblies are produce. Furthermore under the specific case covered in this study, the influence of the velocity of flow seems to not have a significant effect on the ultrasound propagation. This code has been implemented for liquid sodium in the heart of a SFR, but it could be used for many fluids, such as water or oil, and could be used for other kind of ultrasonic instrumentation such as debitmetry or telemetry.

Finally, two experimental verifications of the model used in this code are currently being carried out. In the first one, thermal inhomogeneities are created in a specific fluid and ultrasonic waves are propagated in it. The fluid is chosen for its ability to slow the convection. First results are presented in reference 30. The other one uses Kelvin-Helmholtz instabilities in water with a homogeneous temperature to determine the influence of turbulence on acoustic wave propagation. For these experiments, a comparison between the numerical data from the simulations and the experimental data is realised to determine the validity of the theoretical ray model.