Matthieu Petit
email: matthieu.petit@univ-amu.fr

Lisa Michez

Jean-Manuel Raimundo

Philippe Dumas

Electrical and optical measurements of the bandgap energy of a light-emitting diode

Keywords: numbers: 71.20.Nr, 73.40.Lq, 73.40.-c bandgap, semiconductors, p-n junctions, LED, light-emitting diode

1. Introduction

Context

Semiconductors materials are at the basis of electronics [1]. Most of the electronic devices are made of semiconductors. The functioning of semiconductors devices is completely described by quantum physics in the band theory of electrons in solids. The concept of band structure in the physics of semiconductors often retains an element of mystery for our students of first year of school of engineering (third year of university class) specialized in materials science. Indeed the notions developed in this theory are somewhat abstract. Therefore it is interesting to provide students with simple experiments in order to manipulate these notions. A description of the solid physics of semiconductors materials and of the physics of the semiconductors devices can be found in the books of C Kittel [START_REF] Kittel | Introduction to Solid State Physics john wiley[END_REF] and S M Sze [START_REF] Sze | Physics of semiconductor devices wiley[END_REF]. The band structure of a crystalline solid is the relationship between the energy E of the electrons of the solid and the wave vector k (i.e. momentum). It comes from the Schrödinger equation (eq.1) of an approximate oneelectron problem in a periodic potential U (r) with the periodicity of the crystal lattice:

- 2 2m ∇ 2 + U (r) ψ k (r) = E k ψ k (r) (1)
( r) is a vector of the direct lattice, ψ k (r) the wave function describing the given electron, E k and m its energy and mass respectively. In the case of a periodic crystal, the Bloch's theorem gives the form of the solutions of the Schrödinger equation. There are energy E k domains where no solution exist. These domains are forbidden: electrons cannot have energies with values belonging to these regions. Above and below this energy gap are permitted energy regions called bands. The upper band is the conduction band (CB), the lower band is the valence band (VB) as displayed on fig. 1a). If we called E CB the smallest energy of the conduction band and E V B the highest energy of the valence band, the forbidden energy gap range E g is equal to (eq.2):

E g = E CB -E V B ( 2
)
E g is the bandgap energy, one the most important parameters in semiconductor physics. Each semiconductor material has a specific bandgap. For instance silicon (Si) has a bandgap value of 1.12 electronvolt (eV) and germanium 0.66 eV at room temperature (RT).

External excitations such as temperature, light, electric current can enable the "jump" of an electron from VB to the CB. The energy of this electron increases by a value equal to E g . The "jumping" electron leaves a positive charge called "hole" in the valence band. The excitation creates an electron-hole pair (fig. 1b)). The electron can lose energy and goes down to the VB. The lost energy might be transferred to a photon (radiativ transition) and is also equal to E g as shown in fig. 1c). This phenomenon is the luminescence and is used in optoelectronics devices like light-emitting diode (LED).

Electrical approach for the lab work

Students measure the bandgap of semiconductor materials used in commercially available LEDs. In the case of the lab work and to simplify the problem on an electrical point of view, these devices are described as simple p-n junctions. To get a good understanding of the functioning of a p-n junction, we suggest the reader to read the article of E. Boyes [START_REF] Boyes | [END_REF]. The electrical current I through a p-n junction is given by:

I = I sat e qV k B T -1 (3)
V is the voltage applied to the p-n junction, q the elementary charge, k B the Boltzmann constant and T the temperature. The current-voltage characteristic of an ideal p-n junction is drawn on fig. 2. I sat is a function of T and of the bandgap [START_REF] Sze | Physics of semiconductor devices wiley[END_REF]. It can be written as:

I sat = T 3+ γ 2 e - Eg k B T (4)
γ is a constant depending on a given semiconductor material with a value in the range of 1 to 2. Thus the total expression of the current is:

I = T 3+ γ 2 e - Eg k B T e qV k B T -1 (5)
2. Method and implementation

Equipment

The needed equipments consist firstly in a set of LEDs emitting in the visible light range. LEDs will be generally made of III-V semiconductors. III-V semiconductors are alloys from atoms of the third and fifth columns of the periodic table. They are commonly used in the optoelectronics devices. These devices are powered by a power supply in the range of 0-5 V. A current meter with a resolution of 0.01 µA (measurement uncertainty: ±0.05% + 3digits) and a voltmeter with a resolution of 0.0001 volt (measurement uncertainty: ±0.02% + 2digits) allow to do the measure of the electrical current and applied voltage of the semiconductor devices (fig. 3). A potentiometer with a range from 0 to 10kΩ completes the electrical apparatus require for the lab work. Indeed there are three parameters which can varies in this experimental set-up according to eq.5: the electrical current I, the applied voltage V and the temperature of the electronic device T. To simplify the problem, we choose to fix either the current or the voltage with the variable resistor: adjusting the potentiometer will allow to keep the current or voltage constant. In the following parts, we will work with a constant electrical current. Finally an oven heating up to 100 • C with a temperature regulation and a step control of 1 • C will allow to change the device temperature.

Modelling

In order to calculate E g from the electrical measurements versus temperature, eq.5 must be re-written in a simpler way taking into account the available range of the experimental parameters of the lab work: the applied voltage of the semiconductor device is around 5 V, the temperature is in the range of [RT ; 100 • C]. This modelling part is interesting since students have to adapt their theoretical knowledges about the p-n junction to the specific framework of this pratical work. The logarithm of eq.5 gives the following expression:

ln(I) = 3 + γ 2 ln(T ) - E g k B T + ln e qV k B T -1 (6)
Taking into account the range of the experimental parameters, we can estimate that:

e qV k B T -1 ≈ e qV k B T , (7) 3
+ γ 2 ln(T ) ≈ constant, (8)
ln(I) = constant. (9)
Using eq.7, 8 and 9 in eq.6, we get a linear relationship between the applied voltage and the temperature of the device:

V = C q T + E g q ( 10
)
where C is a constant depending on the chosen fixed value of the current. Thus it is possible to estimate the value of E g by plotting the straight line V versus T (in Kelvin). At T=0 K, V = Eg q is the value of E g directly in eV. Evolution of the applied voltage V versus temperature T for a green LED at two different fixed electrical currents I (30 and 70 µA)

Electrical and optical measurements

To start the electrical measurements, the students have to choose a working value for the electrical current according to the semiconductor device they study. Then the oven is programmed with a temperature ramp. The ramp should be quite low to allow a good thermalization of the electronic component and enough time to take the measures. In fig. 4 we display the measurements coming from four experiments: we have studied a green LED at two different fixed currents (30 and 70 µA). The applied voltages versus the temperature are drawn. As forecast by eq.10, we obtain linear relationships between V and T. Regarding the optical measurements, the light emitted spectrum of the LED is recorded with a USB-650 Red-Tide spectrometer from Ocean Optics [START_REF]Ocean Optics I[END_REF]. This small spectrometer is connected via an USB port to a computer and is very user friendly.

Results and discussion

Electrical measurements versus temperature

The results of the linear fits of our experimental data (fig. 4) are summarized in table 1.

First of all to discuss these results, students should seek what are the semiconductor materials commonly used to manufacture the LED components. Obvioulsy the reference numbers leading to the corresponding device datasheets are relevant to find informations. One of the traditional III-V semiconductor material used to produced green wavelengths is the gallium phosphide (GaP) which has a theoretical bangap energy of 2.26 eV at 300 K. The datasheet of our used LEDs (model L-1503GD from Kingbright) confirmes the nature of the semiconductor [6]. The difference between our experimental values of E g (GaP ) and the theoretical one are 5%. In the next step of the interpretation, the students should provide some explanations to these disparities. Two types of disparities can be contemplated based on experimental and theoretical accounts.

On the experimental point of view, these differences may mainly arise from the temperature measurement: we use the temperature values providing by the oven controller. But these temperatures are overestimated compare to the temperatures of the studied electronic components because of the thermal inertia of the experimental set-up. The influence of this source of errors could be decreased by setting a lower temperature ramp to ensure a better thermalization of the components but the duration of the lab work is limited. On the theoretical point of view the bandgap energy is itself a function of the temperature. This point is not taken into account in our simple electrical model of the p-n junction. Through the Bloch's theorem, E g depends on the lattice parameters of the semiconductors which depend also on the temperature. Varshni propose the following relationship between E g and T [START_REF] Varshni | [END_REF]:

E g (T ) = E 0 - αT 2 T + β (11)
where α and β are constants for a given semiconductors materials. E 0 is the bandgap energy at 0 K. E 0 , α and β are given in table 2 for GaP.

Thus in the temperature range of this practical work, the variation of E g is presented on fig. 5.

The difference of E g between RT and 373 K (100 • C) is around 3%. So this effect may be taken into account to explain the discrepancies between the value of E g found by our electrical measurements versus temperature and the theoretical values given at RT.

T (K)

optical spectrum of the green LED

In the case of the LEDs it is possible to access the (optical) bandgap energy by recording the light spectrum emitted by the device. The spectrum of the used green LED is shown on fig. 6.

The maximum light intensity is recorded for a wavelength λ = (561 ± 5)nm. This corresponds to an energy of E = hc λ = (2.21 ± 0.02)eV which is suitable for a GaP based LED. The reference value of the wavelength given by the manufacturer datasheet is λ ref = 565 nm corresponding to an energy of 2.20 eV.

How to explain the difference between the electrical and optical methods?

The last step of this lab work is to propose an explanation of the difference between electrical and optical values of E g . The main reason is linked to how a LED is really designed. A LED is basically a p-n junction working with a forward bias. Electrons are the majority carriers in the n-type region. They are injected in the p-type region where electrons can recombined with the holes. Thus LEDs are made of semiconductors where the radiative recombinations probability is high. The light colour emitted by the diode depends directly on the bandgap energy of the semiconductor materials where the radiative transition takes place. In order to get a high light emission yield a very good chemical and crystallographic quality of the materials is required. Practically it means that a LED will be built from successive epitaxial (a definition of "epitaxy" can be found on Wikipedia [START_REF]Wikipedia 2015 Epitaxy -wikipedia, the free encyclopedia[END_REF]) stacking layers of semiconductors like shown on fig. 7.

Therefore the LED is not strictly speaking a simple p-n junction that we have used to model the devices for the electrical measurements. The optical measure gives the bandgap energy of the active light emitting layer alone. On the contrary the electrical experiment give access to an "average" bandgap energy of all stacking layers.

Summary

This article provides a quick description of the underlying quantum physics explaining the semiconductor devices. Basic devices such as LEDs can be approximated by p-n junctions. Thus the electrical current flowing through the component depends only on the applied voltage, the temperature and the intrinsic semiconductor bandgap energy E g . A simple electrical set-up combined with an oven allowing the variation of the temperature gives rather good estimation of E g . Moreover the bangap energy of LEDs can also be measure from the wavelength of the emitted light. In this lab work designed for first years university students with a major in materials science and a duration of 3 to 4 hours, we show that by this two methods the E g of GaP used in a green LED is obtained with about a 5% error. The most interesting part is for the students to try to explain theirs results coming from electrical and optical measurements by the various knowledges they have on solid physics and those they find doing some bibliographical researches on the internal design of optoelectronics devices.

Figure 1 .

 1 Figure 1. a) Simplified band diagram of a semiconductor, b) in an excited state with an electron-hole pair, c) when a radiative transition occurs emitting a photon with an energy Eg = hν.

4 Figure 2 .

 42 Figure 2. Ideal current-voltage characteristic of a p-n junction.

Figure 3 .

 3 Figure 3. Electrical scheme of the experimental set-up.

Figure 4 .

 4 Figure 4.Evolution of the applied voltage V versus temperature T for a green LED at two different fixed electrical currents I (30 and 70 µA)

Figure 5 .

 5 Figure 5. Variation of the gap energy with temperature for GaP according to eq.11 and table 2.

λ (nm) 800 Figure 6 .

 8006 Figure 6. Light spectrum of the green LED recorded with a UV-visible spectrometer.

Figure 7 .

 7 Figure 7. scheme of a LED with the electrical contacts in black, the active p-n junction layer and the successive epitaxial layers which guarantee a very good chemical and crystallographic quality of the active layer.

Table 1 .

 1 Results coming from the linear fits of the experimental data (fig.4)

			green LED
			30µA 70µA
	A = C q ×10 -3 ±0.02 × 10 -3	-2.41	-2.28
	B =	Eg q (eV)	2.37	2.37
	±0.06
	R 2		0.998	0.998

Table 2 .

 2 Values of the parameters of eq.11[START_REF] Varshni | [END_REF],[8]

	Semiconductor	E 0 (eV)	α × 10 -4	β
	GaP	2.34	6	460


John Brandon Lowe

Computerized tools for reconstruction in Tibeto-Burman

Proceedings of the Fifteenth Annual Meeting of the Berkeley Linguistics Society, Feb 18-20 1989: 367-378


As it is easy to create scenarios containing various static and moving obstacles, evaluation of control algorithms in several situations is its main feature. The adopted architecture for the simulation tool also enables evaluation of road traffic scenarios constructed on various levels of interaction or cooperation among cybercars.

Index Terms-Cybernetic Transportation System, Driverless

Vehicles, Cybercars, Simulation.

I. INTRODUCTION

he common use of personal vehicles all over the world causes safety, congestion and pollution problems, especially in large cities. To solve these problems, an innovative approach, called Cybernetic Transportation System (CTS), has been proposed [START_REF] Parent | Automated public vehicles: A first step towards the automated highway[END_REF][START_REF] Parent | Cybercars: Past, Present and Future of the Technology[END_REF]. CTS is a city transportation system based on a fleet of unmanned electrically powered carlike robots. These cybercars, equipped with sensors, have the ability of autonomously transporting people in an urban environment. The highly dynamic nature and unpredictability of such environments makes the development of a cybernetic transport system very complex. Therefore, the use of a simulator appears to be a good solution to facilitate developments of CTS driving capabilities as it allows multiple testing of the same road traffic scenarios and easy evaluation of the overall system's traffic performance.

This work was supported by European Commission Information Society

Technologies under Grant 028262 and Griffith University.

Sébastien Boissé is with IMARA Team, INRIA Rocquencourt, France and with Ecole Centrale Nantes, Nantes, France (e-mail: Sebastien.Boisse@eleves.ec-nantes.fr).

Rodrigo Benenson, Laurent Bouraoui and Michel Parent are with IMARA Team, INRIA Rocquencourt, France (e-mail: Rodrigo.Benenson@inria.fr, Laurent.Bouraoui@inria.fr, Michel.Parent@inria.fr).

Ljubo Vlacic is with ICSL Laboratory, Griffith University, Brisbane, Australia (e-mail: L.Vlacic@griffith.edu.au) Lots of simulation software packages have been developed in the past and more recently, some 3D simulators were released, again reducing the gap between simulation and reality. Among all of those simulators, some of them deal with wheeled robots, including various sensors like cameras, range sensors or contact sensors. Most of these tools are built upon Open Dynamic Engine (ODE) [4], an open source library providing high performance for rigid body dynamics simulation.

In this area, two free tools, OpenSim [1] and Gazebo [START_REF] Koenig | Design and Use Paradigms for Gazebo, An Open-Source Multi-Robot Simulator[END_REF] are able to simulate several wheeled vehicles with sensors in a 3 dimensional environment. The OpenSim provides real-time rendering for simple shaped robots while Gazebo can mimic a complex world thanks to simulated vehicle sensors. The Webots [START_REF] Michel | WebotsTM: Professional Mobile Robot Simulation[END_REF], a commercial simulation software for mobile robots proposed by Cyberbotics, simulates several models of vehicles and various sensors. This prototyping application permits testing of control algorithms on various simulated robots.

In another area, USARSim [START_REF] Koes | Validating USARsim for use in HRI Research[END_REF] has been developed to be the basis for the RoboCup Urban Search And Rescue (USAR) simulation competition. This simulator, based on the Unreal Tournament Game Engine, has the particularity to simulate a camera with accuracy thanks to its high-quality rendering.

Finally, CycabTk [2] has been developed at INRIALPES (France) to simulate a cybercar with its sensors (camera, SICK laser, …), using the same interface as real vehicles. This simulator runs under linux and uses mgEngine as 3D engine [5].

While some simulation tools tend to be as generic as possible in order to simulate several robots and sensors, we have chosen to develop a simulation tool that can accurately mirror the vehicle we are using at INRIA (Lara Research unit, France) and Griffith University (Intelligent Control System Laboratory, Brisbane, Australia), in a 3D virtual environment. In addition, to appropriately simulate CyCabs (name of the cybercars developed by Robosoft) dynamics and sensors, this tool is designed to face CTS issues by providing an opportunity to test the behavior of a fleet of vehicles in the same simulated environment. The latter is composed of a world model (map containing roads, buildings and other static obstacles), moving obstacles and cybercars. The 3D models of both maps and moving obstacles are fully customizable as T they are created with 3D modelers. In this way, we will be able to deal with interactions between vehicles involved in many different traffic situations.

In section II, we present the architecture of the framework, and then give details about the simulation in sections III and IV. In section V we present the different traffic situations considered. The first results obtained with this simulator are highlighted in section VI and finally, limitations and future works are discussed in section VII and section VIII.

II. FRAMEWORK ARCHITECTURE

Due to the high dynamicity of urban environments, the cybercar control software program needs to obtain data from sensors and to analyze them at a high frequency, which is computationally expensive. This often requires the whole CPU of a powerful computer to work properly. Moreover, each cybercar is controlled by its own computer; therefore, we adopted the same architecture when designing the simulation tool. The simulator acts as a server and each control software program (clients), which can be run on distinct computer, is connected to a simulated vehicle inside the simulator through a TCP/IP socket (Figure 1). Basically, the simulator sends the simulated sensors' data to each control software program and receives in return commands from them to make the vehicle move inside the virtual environment (Figure 2).

The design of this framework allows the different remote control software programs to communicate between each other through the network without transmitting by the simulator, using either Ethernet or Wifi links. In this way, the control software program, used to exchange data between cybercars, can also be tested with this simulator.

Remote control software programs can be easily connected to the simulator thanks to Zero Configuration Networking (ZeroConf) standards [START_REF] Guttman | Autoconfiguration for IP Networking: "Enabling Local Communication[END_REF], permitting automatic detection of the server's address (IP + PORT) over the network (Figure 3).

The simulator, running under windows, uses a crossplatform interface allowing control software program developed under windows or linux to connect to it. The interfaces to access sensors' information and to send commands are the same as in real vehicles, so the same application can be run on both real and simulated CyCabs.

A. Simulation loop

To simulate the dynamics of the environment, we use a physics engine able to compute interactions between all objects introduced into the virtual world. Open Dynamic Engine (ODE), a free library for simulating rigid body dynamics allows us to introduce vehicles easily in a 3D environment.

Each simulator requires an interface to show the results of its computations. Rendering of 3D objects appears to be the best solution since we wish to see the results of the simulations in real-time, and in a format than can be easily interpreted by people with non-scientific backgrounds. Further, it makes it easy to save the results in a video file to display our results during exhibitions.

However, rendering is just a way to see the results and must not affect the simulation itself. To render 3D textured objects is not an easy task; it sometimes requires significant CPU time and some care must be taken to minimize the impact on other tasks.

To fulfill these requirements, the simulator has been built using the 3impact Engine [3], which is a low cost library integrating an ODE-based physics engine and a rendering engine. We made this choice because the rendering engine operates without affecting the physics engine. The frequency of the physics engine can be set (75Hz is used as it provides sufficient accuracy for computations), and will remain constant during the simulation. The simulation main loop, computing all the motions and checking for collisions, will always be executed at the same frequency and the rendering will be done only if time is available so the frame rate will decrease if the computer is not fast enough.

B. Models

The 3impact Engine does not require any specific 3D model editor; any modeler able to export to DirectX format can be used (for example 3D Studio Max [15], LightWave 3D [16], MilkShape 3D [17], AC3D [18]). Therefore there are no limitations on the physical structure of the models loaded into the simulator.

1) Vehicle models

The simulated vehicle (CyCab) is defined by three 3D models:

a) The first model, made only with triangles, defines the shape of the vehicle. It also contains texture-mapping information. Created from the CAD model, it is only used for rendering so the number of triangles has been reduced (around 7000) to avoid losing too much time when rendering it (Figure 4a).

b) The second model, made only with sphere primitives is the collision model. All collision tests between the vehicle and something else will be performed using this sphere-based model. Therefore, this model must match closely the first model. This model is built up with 41 spheres (Figure 4b).

c) The third model, a simple box bounding the car, defines the mass distribution of the vehicle.

The combination of the last two models defines the model used by the physics engine.

Once the models of the vehicle have been introduced into the virtual world, four independent wheels are joined to this vehicle with strength and damping suspension parameters. Friction parameters between the wheels and the ground are also defined, allowing us to simulate realistic performance of the vehicles [3].

2) Moving obstacle models

For moving obstacles, we use the 3 models described above to define a vehicle. Collision checks are performed using the sphere-based model, so the accuracy of the simulation is dependant on how this model is close to the mesh model (number and size of spheres).

3) World model

World models are 3D maps representing a real area and are only composed of static obstacles (roads, buildings, parked cars…). These triangle-based models, containing information about textures, are used to compute collision checking. However, in order to avoid time-expensive collision tests, world models should only contain basic shapes made with few triangles. Because collision checks between rays and triangles are slower than collision checks between rays and spheres, more detailed static models should be introduced as moving obstacles with null speed in order to use a sphere-based model for collision tests.

4) The whole virtual environment

The various objects (vehicles, world model, moving obstacles) are introduced separately into the simulator. To process the simulation of the global environment, interactions between those objects have to be defined. These interactions can be defined with collision couples specifying friction parameters between two different objects. Creating collision couples between all objects is not necessary and would necessitate a lot of time; so only couples between vehicles and the map are added (only CyCabs need accurate simulation). There are five couples per vehicle (one for the vehicle body, and one for each wheel).

C. Rendering

A camera is introduced into the virtual environment to get a view of the simulation; the position of the camera can be modified in runtime with the keyboard. We can focus on each of the different simulated vehicles from various directions and also are able to get a global view of the scene.

The DirectX-based rendering engine renders the different object seen by the camera in the 3D environment.

IV. SIMULATION OF SENSORS AND ACTUATORS

A. Sensors

To get a good simulation, models of the sensors we are using have to be as accurate as possible. Moreover, we need to use the same interface to access the virtual sensors' data and the real sensors; we can then use the same application to drive the real vehicle and the simulated vehicle.

The time available to analyze sensors' data and compute commands is dependent on the frequency with which sensors send data, because for each step in the calculation, we want to take new data into consideration. Therefore, it is essential that the simulator provides data at the same frequency as in reality.

As a starting point, three sensors have been implemented: Laser, GPS and compass.

1) Laser sensor

The laser we are simulating (IBEO) scans the world in 181 directions, each separated by one degree. The laser sensor returns a vector of distances obstacles, each element corresponding to a given direction and we check whether there is any intersection between the rays and obstacles. The laser's range has been set to mirror the real laser.

The different types of obstacles are treated separately to improve calculation speed. In the first step, moving obstacles and simulated vehicles are considered as single bounding spheres; if a collision is detected between a ray and one of these spheres, the accurate sphere collision model of the corresponding object is used to determine the accurate distance to the obstacle. For each direction, if more than one obstacle is detected, the shortest distance is returned.

By default, the scanning frequency of the laser is 8 Hz but this value can be adjusted.

2) GPS sensor

To simulate the GPS sensors, we approximate the small area of the globe where we are located by a map with a 2D system of coordinates. By default, the reference point of this system is the geometric centre of the virtual world model, but we can translate and rotate this coordinate system to match a virtual map with a real environment.

GPS simulated sensors return (X, Y) coordinates relative to this coordinate system. The default frequency of this sensor is 1 Hz but is also adjustable.

3) Compass sensor

Compass sensors give the orientation of the vehicle, relative to north as defined in the 3D model of the world by the simulator.

4) Noise

To make our modeling more realistic, some noise (with customizable Gaussian distribution) can be added to the sensors' data sent to the client.

B. Actuators

On the real CyCab, each wheel is controlled independently by an electric motor. It is the role of the control software program to determine feasible commands before sending them to the actuators as the CyCab has limits on acceleration and wheel angle. In the simulator, the command (speed in m/s and steering angle in radians) will also be applied to each wheel while keeping in mind the physical constraints of the real CyCab. If a non-feasible command is sent, the highest possible command will be applied.

Each command sent by the client to the simulator is temporarily stored and applied to the intended CyCab until new commands are received by the simulator. One of the most important features of this software is that it can easily model various traffic scenarios. For instance, in addition to the static environment, it is possible to add moving obstacles like a pedestrian crossing the road or walking on the side of the road.

The trajectory of moving obstacles can be customized with a list of waypoint to reach the map coordinates system and the speed. These obstacles are moved in runtime by the simulator and do not require any external software to pilot them. A simple control system inside the simulator dictates the orientation of the moving obstacle towards its next waypoint. The simulator applies the specified speed and directs the object to its next valid waypoint.

A small library of maps has been modeled, ranging from very simple cases to more complex maps containing buildings and other static obstacles. They typically represent various traffic conditions that may be encountered in reality like different types of roads (straight or curved) or intersections ("T", "x"). With all of these models, many scenarios can be set by customizing:

1) The number of simulated vehicles 2) Their initial location 3) The number of moving obstacles 4) The type of moving obstacles 5) The predefined motion of these obstacles Each scenario is defined in a formatted file where all the parameters described above are stored. Figure 5 schematizes the simulation of the virtual environment and Figure 6 shows the view of a scenario involving several simulated vehicles and pedestrians in a city-like crossing situation.

During the simulation, virtual vehicles are driven by remote applications and execute maneuvers defined by control algorithms, taking into consideration the different obstacles detected by the sensors in the virtual world. The camera introduced inside the simulator allows us to analyze the cybercars behavior. The real and simulated CyCabs are not completely isolated from one another and can interact with each other. The level of this interaction depends on the capabilities of the control software program to exchange data. This interaction is possible because the simulation runs in real-time (provides sensors information at the same frequency as real CyCabs), processes data in the same way as the real CyCabs and the world model can be created to mirror the area around the real CyCabs, including GPS location and orientation of the CyCabs. In this way, simulations involving virtual and real vehicles can be run before attempting the same scenario with only real vehicle.

VI. RESULTS

The simulator is relatively new but has already proved its capability to test new control algorithms by simulating CyCabs' driving performance.

The simulator was originally developed to simulate the interaction of cybercars at crossings. The goal was to evaluate whether an exchange of data between two vehicles at a crossing would improve their ability to avoid each other by anticipating the behavior of the other.

Based on a control software program [START_REF] Benenson | Integrating Perception and Planning for Autonomous Navigation of Urban Vehicles[END_REF] able to drive the car to its final destination while avoiding static and moving obstacles [START_REF] Petti | Safe motion planning in dynamic environments[END_REF], the new control software program now contains communication capabilities to broadcast data through a UDP socket over a wireless network. This communication module allows sending information about various obstacles that have been perceived by the sensors and planned trajectory of the vehicle. This module also receives the broadcasted UDP packets and adds them to the world model previously created with its own perceived obstacles (Figure 7).

The vehicle's world model now contains the obstacles perceived by the sensors, the obstacles perceived by the other cybercars and their last computed trajectories represented as a succession of obstacles in time. Therefore, the planning module will consider all of these obstacles to compute its new trajectory. The simulator was first used to make the existing algorithm from the real CyCab more robust and stable. The first tests of the new algorithm were run in a scenario called "face to face encounter" that involves two vehicles facing each other set on a collision course.

The results of these tests show that the vehicles anticipate and avoid the collision by adapting their trajectory to the new received future obstacles (i.e. the trajectory of the other cybercar). Without the communication, tests had shown that the vehicles can stop in front of each other, unable to move anymore. Then, without any priority rules in the planning algorithm, vehicles were able to pass beside each other. Now that the algorithm has been validated, two real CyCabs are being equipped with the appropriate sensors to run the same experiments for real.

VII. LIMITATIONS

The purpose of this simulation software is to test functionalities of autonomous driving algorithms in city-like road traffic situations, helping design and development of Cybernetic Transportation System. However, the presented simulation tool cannot be used for the purpose of evaluating the global performance of CTS systems. Another tool will be needed to accurately monitor the traffic (average time to reach a specified location, traffic flow in crossings, etc.) and prove that CTS can resolve congestion problems in cities.

To remove unessential calculations, only collision couples between vehicles and the static world have been defined. This allows more simulated vehicles and more moving obstacles to be added. Therefore, it is not possible to get accurate feedback about collisions between other objects, and the results of the simulation (obstacle avoidance for instance) have to be visually interpreted through the rendering window.

As previously mentioned, the main loop frequency remains constant during the simulation ensuring the accuracy of the simulation. Therefore, the number of CyCabs able to be introduced to any given scenario is heavily dependant on the processing power of the computer running the simulator. The quick collision test capabilities of 3impact engine (thanks to sphere-based models) let us introduce up to six simulated vehicles in a world up to 500*500 m2 containing static and moving obstacles on a Pentium 4 2.8 GHz with 512 MB of RAM. If the computer running the simulator is not fast enough to process the main loop at the set frequency, a black screen is displayed (rendering frequency tends towards zero). Therefore, a more powerful computer is needed to simulate more complex scenarios.

VIII. CONCLUSION AND FUTURE WORK

The simulation software presented in this paper aims to help development of CTS by simulating several cybercars with physics and sensors in a virtual 3D environment. 3D models of worlds and moving obstacles are fully customizable as most 3D modelers can be used to create them. Therefore, it is possible to build very specific environments (traffic situation, crossings, pedestrian…).

The real vehicles are equipped with more sensors (camera, odometers, …) than are currently simulated. Therefore, to keep the simulation software up-to-date, further development will include the addition of new sensor simulation models.

In the meantime, in order to simulate a complete traffic scenario including communication between cybercars, some features may be added to make communication more realistic. The current simulator setup permits communication between all vehicles being simulated, regardless of whether they would be able to communicate in real life (maximum wifi range), as they are all connected via ethernet or WiFi. New modules will be added to emulate the WiFi structure regarding of the relative location of the vehicles in the virtual world and to filter packets that actually can not reach their target.

Moreover, first experiments have shown that cybercars which use the Optimize Link State Routing (OLSR) protocol [START_REF] Bouraoui | OLSR improvement for distributed traffic applications[END_REF] allow for an ad-hoc network to be created between the vehicles, which dynamically updates itself with respect to the relative position of the vehicles and the radio links present between them. As communication between cybercars seems to be a very important research topic for the next few years in CTS development, a simulation of this protocol based on information available in the simulator could be interesting as well.

Cybernetic

  Transportation Systems Design and Development: Simulation Software Sébastien Boissé, Rodrigo Benenson, Laurent Bouraoui, Michel Parent, Ljubo Vlacic Abstract-The growing number of vehicles saturates cities in terms of congestion and pollution. A Cybernetic Transportation System (CTS) appears to be a way to resolve those problems. Based on a network of clean driverless vehicles (the cybercars), CTS aims to improve safety and organization of urban transport by providing a door-to-door complement to efficient and fast mass transport. The simulation software presented in this paper has a goal to facilitate the development of such a transportation system. It simulates several cybercars in a dynamic virtual 3D environment, and provides sensors information in real vehicles.

Figure 1 :

 1 Figure 1: Framework architecture

Figure

  Figure 2: Simulation

Figure 3 :

 3 Figure 3: Zero Configuration Networking

  Figure 4: CyCab models

Figure 5 :

 5 Figure 5: Virtual environment simulation

Figure 6 :

 6 Figure 6: scenario example

Figure 7 :

 7 Figure 7: Remote control software program architecture


Bassem Barhoumi
email: bassembarhoumi@hotmail.fr

Boundary Element Method Formulation for axisymmetric acoustic problems in a subsonic uniform flow

Keywords: Convected Helmholtz equation, modal Green's function, axisymmetric boundary element method, singular integrals, axisymmetric duct

This paper presents a new analysis method and numerical development of the direct boundary element method formulation for axisymmetric acoustic radiation and propagation problems in a subsonic uniform flow. It governed by the homogeneous axisymmetric convected Helmholtz equation independent of the flow direction based on the Fourier coefficient, the three-dimensional convected Green's function in free space according to the azimuth angle represented by the modal kernel arising in the acoustic sources. This axisymmetric integral formula is expressed only in two new operators; the first, concerned the particular normal derivative similar to the temporal derivative and the other, concerned the convected normal derivative reduce the convective terms of normal derivative and the flow direction derivative used in the classical formulations developed in the literature. When, the main source is taken on the generator, the integration modal Green kernel and its convected derivative are performed partly analytically in terms of Laplace coefficients evaluate by a simple recursive formula and partly numerically using a Gauss-Legendre quadrature standard formula to facilitate the analytical resolution and the conventional numerical problems. The boundary element method formulation used to illustrated and tested by applying the infinite axisymmetric cylindrical duct in a subsonic uniform flow.

Introduction

Noise nuisance caused by airplanes, satellites and automobiles are becoming a major and growing concern. They are generated by the radiation or the propagation of acoustic waves in a fluid medium in flow or at rest. Thus, several methods have been developed to help understand the physical phenomena involved to reduce this noise such as, the finite element method and boundary element method, when these domains are in revolution, it is necessary to develop new representations for these methods. Firstly, The FEM in axisymmetric problems (FEMA) is adapted in bounded interior region where [START_REF] Craggs | A Finite Element Method for Damped Acoustic Systems: An Application to Evaluate the Performance of Reactive Mufflers[END_REF][START_REF] Zhang | A coupled FEM/BEM formulation for acoustic radiation in a subsonic non-uniform flow[END_REF] are determined the acoustic radiation and propagation respectively, for spherically and cylindrically symmetric problems in uniform flow. This simplification by Prandtl-Glauert transformed the flow axisymmetric domain to a no-flow case [START_REF] Balin | Boundary Element and Finite Element Coupling for Aeroacoustics Simulations[END_REF]. However, this method is poorly adapted to resolve the acoustic problems in an unbounded axisymmetric exterior domain. Since, it was been induced reflections or radiation condition at infinity [START_REF] Sommerfeld | Partial differential equations in physics[END_REF][START_REF] Mosher | the influence of a wind tunnel on helicopter rotational noise[END_REF]. This case is better to use the technique of the axisymmetric boundary element method [START_REF] Juhl | An axisymmetric boundary element formulation of sound wave propagation in fluids including viscous and thermal losses[END_REF][START_REF] Beldi | A simple and effective axisymmetric convected Helmholtz integral equation[END_REF] derived from the three-dimensional integral formulation who can be formulated in the transformed acoustic medium based on the Prandtl-Glauert transformation [START_REF] Balin | Boundary Element and Finite Element Coupling for Aeroacoustics Simulations[END_REF] or in the original flow acoustic domain [START_REF] Beldi | Some new results for the study of acoustic radiation within a uniform subsonic flow using boundary integral method[END_REF][START_REF] Wu | A direct boundary integral formulation for acoustic radiation in a subsonic uniform flow[END_REF] and in the no-flow domain [START_REF] Kuijpers | An improved acoustic Fourier boundary element method formulation using fast Fourier transform integration[END_REF].

The main disadvantages of conventional axisymmetric integral formula are caused by the complexity of boundary conditions transformation to the no-flow case, the convection effects of the normal derivative and the flow direction derivative [START_REF] Beldi | Some new results for the study of acoustic radiation within a uniform subsonic flow using boundary integral method[END_REF]. Also, singularity problems of modal Green function and its derivatives can be evaluated by an analytical method [START_REF] Wu | A direct boundary integral formulation for acoustic radiation in a subsonic uniform flow[END_REF][START_REF] Seybert | A special integral equation formulation for acoustic radiation and scattering for axisymmetric bodies and boundary conditions[END_REF] and by fast Fourier transform techniques [START_REF] Kuijpers | An improved acoustic Fourier boundary element method formulation using fast Fourier transform integration[END_REF] containing the additional terms.

In this paper i present a new method to obtain the conventional Boundary integral formulas for axisymmetric acoustic medium (BEMA) based on the axisymmetric convected wave equation and its fundamental solution such as convected modal Greens function. This conventional formula converted into a simple form similar to the no-flow case, dependent on two new operators physically reduce the convective terms and facilitates the analytical resolution of singularity problems. The numerical implementation of BEMA is developed for acoustic propagation in an axisymmetric cylindrical duct and validated by comparison with FEMA method and the analytical solution.

Axisymmetric formulation

The axisymmetric physical acoustic medium can be divided into three-regions; the axisymmetric interior and exterior domain or the propagation and radiation problem in a subsonic uniform flow and the interface domain in a subsonic non-uniform low-flow similar to the regions of an aero-engine contains the source region, the duct domain and the Farfield region.

We neglect the non-uniform flow effect, this acoustic domain transform into two regions Ω , with = i for the interior domain and = e for the exterior domain R 2 in a compressible isentropic fluid, inviscid, of the density ρ ∞ , of the speed sound c ∞ and the flow Mach vector M ∞ = M ∞ e z in the z-direction of a cylindrical coordinate system (r, θ, z) where θ is the revolution angle [START_REF] Sugimoto | Modelling of flow effects on propagation and radiation from bypass ducts[END_REF] are limited by the fictitious closed generator without needed to impose boundary conditions and unchanged the direction flow (see Figure 1). where the signs (+, -) designates the acoustic exterior and interior domain in a subsonic uniform flow and ∂Ω = Γ ∞ ∪ Γ is the generator limiting the axisymmetric domain Ω with Γ ∞ the exterior artificial boundary of the circular generator for sufficiently large radius and Γ is the interior artificial boundary. Which are closed generators of the outward unit normal vector n q (n rq , n zq ) and of the unit tangent vector t q (t rq , t zq ) at the boundary point q.

Axisymmetric wave equation

A linear combination of fluid mechanics equations based on the mass conservation, the linear momentum and the adiabatic state to produce the standard form of the wave equation, to describe the three-dimensional homogeneous convected Helmholtz equation for the harmonic acoustic pressure perturbation is expressed in the form P = p exp(+iω t) of angular frequency ω and p represents the complex amplitude, i is the imaginary unit, which can be written as

H 3D (p) = ∆ 3D (p) + k 2 ε p -M ∞ε .∇ 3D 2ik ε + (M ∞ε .∇ 3D ) (p) = 0, inV ε (1)
where V is the three-dimensional acoustics problems, H 3D signifies the three-dimensional convected Helmholtz operator, ∆ 3D and ∇ 3D are the operators of Laplace and of Gradient in three-dimensional, k = ω /c ∞ is the original wavenumber. The axisymmetric acoustics problems are handled by expanding the physical quantity into a complex Fourier series according to the revolution angle

p (r, θ, z) = n=+∞ n=-∞ p n (r, z) e inθ (2)
where p n is the modal acoustic pressure or the Fourier coefficient of the three-dimensional acoustic pressure of azimuthal order n independent on the revolution angle but still dependent on r and z can be represented in the form

p n (r, z) = 1 2π 2π 0 p (r, θ, z) e -inθ dθ (3)
By expanding all variables of the wave equation Eq. ( 1) in complex Fourier series Eq.

(2) and integrate the resulting equation multiply by exp(-imθ)/2π over the revolution interval [0, 2π] with m ∈ Z, the three-dimensional convected Helmholtz operator H 3D transformed to the axisymmetric convected operator H axi for modal acoustic pressure decompose into no-flow part and a convection part

H axi (p n ) = ∆ + k 2 ε -(n 2 /r 2 ) (p n ) -(M ∞ε .∇) (2ik ε + (M ∞ε .∇)) (p n ) = 0 ( 4
)
where ∆ and ∇ are the operators of Laplace and of Gradient in axisymmetric problems r-z. The fundamental solution of this axisymmetric wave equation for an acoustic source in free space is the convected modal Green's function G k n (m, q), the Fourier coefficient of the three-dimensional convected Green function G k (M, Q) in a unbounded medium at the locations observer-source points are respectively denoted by M(r M , θ M , z M ) and Q(r Q , θ Q , z Q ) identified in axisymmetric problems by this points m(r m = r M , z m = z M ) and q(r q = r Q , z q = z Q ) (Figure 2).

H axi (G k n (m, q)) = ∆ q + k 2 ε -(n 2 /r 2 q ) G k n (m, q) +(M ∞ε .∇ q ) 2ik ε -(M ∞ε .∇ q ) G k n (m, q) = - 1 2πr q δ(q -m) (5)
where ∇ q is the axisymmetric gradient operator at point q and δ is the Dirac delta function. The intrinsic form of three-dimensional convected Green function in free space G k (M, Q) for M = Q independent of the flow direction, which is represented by

G k (M, Q) = exp -ik * MQ.M * ∞ + R * /(4πcR * ) ( 6
)
where c is the Prandtl-Glauert factor, k * = k /c is the convected wavenumber, M * ∞ = M ∞ /c is the convected Mach vector and R * is the three-dimensional convected radius, given by

R * = R 2 + (MQ.M * ∞ ) 2 (7)
where R = MQ is the physical distance between source-observer in three-dimensional.

Using the cylindrical coordinate system of the points M and Q, the convected radius R * can be rewritten as a function of the azimuth angle

R * = τ * 1 -τ 2 cos 2 (θ) (8)
where the angle θ = β/2 and the parameters τ and τ * given by

τ 2 τ * 2 = 4r q r m , τ * = r * 2 + 4r q r m (9)
The convected distance r * = r 2 + (mq.M * ) 2 is the axisymmetric convected physical radius such as r = mq is the axisymmetric physical distance between axisymmetric observersource points m and q.

Boundary conditions

For the physical acoustic interior and exterior problems it is necessary to introduce the boundary conditions of the fictitious closed generator and the exterior boundary at infinity represented by the duct conditions and the Sommerfeld radiation condition in a subsonic uniform flow. Firstly, the duct conditions can be represented by the rigid duct or vibration duct without the need to apply. Secondly, the Sommerfeld radiation condition represents the outgoing radiation waves, can be rewritten the radiation condition in a subsonic uniform flow [START_REF] Balin | Boundary Element and Finite Element Coupling for Aeroacoustics Simulations[END_REF][START_REF] Mosher | the influence of a wind tunnel on helicopter rotational noise[END_REF][START_REF] Beldi | Some new results for the study of acoustic radiation within a uniform subsonic flow using boundary integral method[END_REF] take the advanced form, given by lim R→∞ R dp dR + ik ∞n p)} = 0 [START_REF] Kuijpers | An improved acoustic Fourier boundary element method formulation using fast Fourier transform integration[END_REF] where d/dR is the particular temporal derivative in harmonic-time, k ∞n = k t -kM ∞ n is the non-standard convected wavenumber depended the mean flow with k t is the tangential convected wavenumber transformed by Prandtl-Glauert of the tangential Mach number

M ∞ t = M ∞ .t q , and M ∞ n = M ∞ .n q is the normal Mach number, given by kt = k ε 1 -M 2 ∞εt ( 11
)
This advanced form of the Sommerfeld radiation similar to the radiation condition for a fluid at rest without transforming the domain and the non-standard convected wavenumber k ∞n becomes the original wavenumber k , the particular temporal derivative d/dR becomes the normal derivative ∂/∂R

lim R→∞ R ∂p ∂R + ik ε p)} = 0 (12)
These two Sommerfeld radiations conditions are verified by the convected Greens function in free space.

Boundary element method for axisymmetric problems

The Boundary element method for axisymmetric acoustics problems is obtained by multiplying the homogeneous axisymmetric convected Helmholtz equation Eq. ( 4) for the modal acoustic pressure H axi (p n (q)) at an point q by the convected modal Green's function G k n (m, q) and multiplying the inhomogeneous convected Helmholtz equation Eq. ( 5) for the convected modal Green's function H axi (G k n (m, q)) by the modal acoustic pressure p n (q), subtracting these two terms followed by integration over the axisymmetric domain Ω with respect to q = m can be written as

Ωε p n (q) 2πr q δ (q -m) dΩ ε q = Ωε ∆ q G k n (m, q) p n (q) -∆ q p n (q) G k n (m, q) dΩ ε q -2ik ε Ωε G k n (m, q) (M ∞ .∇ q ) p n (q) + p n (q) (M ∞ .∇ q ) G k n (m, q) dΩ ε q - Ωε G k n (m, q)M ∞ .∇ q (M ∞ .∇ q )p n (q) -p n (q)M ∞ .∇ q (M ∞ .∇ q )G k n (m, q) dΩ ε q (13
) Using the Dirac property for the left term, the axisymmetric Greens identity formula for the right first term, applying the Gradient property on the right second term and the divergence theorem for the last term. The surface integral for the right terms in Eq. ( 13) transformed into a boundary integral for the modal acoustic pressure can be expressed as

p n (m) 2π = ± ∂Ωε p n (q) ∂ nq G k n (m, q) -G k n (m, q) ∂ nq p n (q) r q dΓ q ± ∂Ωε 2ik ε M ∞εn p n (q) G k n (m, q)r q dΓ q ± ∂Ωε M ∞εn G k n (m, q) (M ∞ε .∇ q ) p n (q) -p n (q) (M ∞ε .∇ q ) G k n (m, q) r q dΓ q (14
) where ∇ q = σ n n q + σ t t q is the axisymmetric Gradient operator such as the couple (σ n = ∂ nq , σ t = ∂ tq ) designates the normal and tangential derivative with respect to the direction indicated by the vectors n q and t q , respectively.

In a no-flow case, the axisymmetric boundary integral equation represented by the right first term of Eq. ( 14) depending on the symmetric modal Green function of G k = exp(ikR)/4πcR and its antisymmetric modal normal derivative ∂ n Q (G k ), but in a uniform flow case of the integral formula Eq. ( 14) they becomes more complicated because their contains of convection terms on the right-hand side of the normal derivative ∂ nq and the flow direction derivative (M ∞ .∇ q ) for the convected modal Green's function. These convection operators can be converted into a simple particular normal derivative d nq take into account of the flow direction similar to the particular temporal derivative adapted to the theoretical developments and numerical implementation, defined by

d nq (•) = ∂ nq (•) -M ∞εn M ∞ .∇ q (•) (15)
The particular normal derivative of the convected modal Green function

d nq (G k n ) is the Fourier coefficient of d n Q (G k ) = ∂ n Q (G k ) -M ∞εn (M ∞ε .∇ Q )(G k
) with respect to the direction indicated by the unit outward normal vector n Q at the boundary point Q, of order n according to the azimuth angle β. Using the normal derivative of the convected Greens function ∂ n Q (G k ) and the flow direction derivative(M ∞ .∇ Q )(G k ) given by [START_REF] Beldi | A simple and effective axisymmetric convected Helmholtz integral equation[END_REF], one obtains

d n Q G k (M, Q) = -(1 + ik * R * ) r nβ R * 2 + ik ε M ∞εn G k (M, Q) (16)
where r nβ = MQ axi .n q is the normal distance depended on the azimuth angle , the vector

MQ axi is coordinate (r Q -r M cos(β), z Q -z M ).
The particular normal derivative of the convected Greens function Eq. ( 16) decomposed into two terms, the first way, it generalizes the normal derivative expression in no-flow case and the second way it explains by the convection influence. Substitution of the particular normal derivative Eqs. ( 15) and ( 16) into the axisymmetric integral formula Eq. ( 10), yields a relatively simple form of the modal acoustic pressure

p n (m) 2π = ± ∂Ω G k n (m, q) d nq p n (q) -p n (q) d k nq G k n (m, q) r q dΓ q (17)
These modified axisymmetric integral formula Eq. ( 17) form an explicit expression for the modal acoustic pressure and its particular normal derivative becomes as a no-flow integral formula given by the right first term of Eq. ( 14) over the generator used to determine the modal acoustic pressure in the exterior and interior domain. From the particular normal derivative operator Eq. ( 16), we obtain the convected normal derivative of the convected modal Green function

d k nq G k n = d nq G k n +2iM ∞εn G k n as the Fourier coefficient of d k n Q (G k
) of order n and the azimuth angle β, given by

d k n Q G k (M, Q) = -(1 + ik * R * ) r nβ R * 2 -ik ε M ∞εn G k (M, Q) (18)
The axisymmetric integral formula expressed by Eq. ( 17) over the large circular generator Γ ∞ represent the contribution wave on the interior due of the source placed at infinity such as q → ∞ implies that the integration is properly convergent and becoming zero value. This integral term leads two convected boundary conditions verified by the modal acoustic pressure and the modal Greens function. The first indicates that the pressure must decrease, as we go to infinity (non-reflection condition) and the second generalize the usual radiation condition in a subsonic uniform flow, ensures that all radiated waves are outgoing justified by the choice of these new operators Eqs. ( 15) and [START_REF]Acoustics of ducts and muffler[END_REF]. When the point m → q, the modal kernel function G k n (m, q) and their convected normal derivative contain singular part localized as well as regular part more global.

Singular behavior

A classical procedure to separate these singular and regular parts obtained by a Taylor expansion, decomposing the convected Greens function G k = G 0 +g k into a static singular part G 0 contain a singularity at M → Q and regular part g k depending on the wavenumber

G 0 (M, Q) = 1 4πcR * (19) g k (M, Q) = exp(-ik * R) -1 4πcR * (20)
where R = MQ.M * ∞ + R * is the special radius. The convected normal derivative from this function to the same form T k = T si k + T re k where the singular function T s k i = d k n Q (G 0 ) is convected normal derivative of the static Green function given by Eq. ( 6) and the function T re k is the regular part, given by

T si k (M, Q) = - r nβ 4πcR * 3 + ik ε M ∞εn G 0 (M, Q) (21) T re k (M, Q) = 1 -(1 + jk * R * ) exp(-ik * R) r nβ 4πcR * 3 + ik ε M ∞εn g k (M, Q) (22)
In the same way, the Fourier coefficients of the convected Greens function and its convected normal derivative in the integral equation Eq. ( 17) decomposed as follows

G k n (m, q) = G 0 n (m, q) + g k n (m, q) (23)
T k n (m, q) = T si kn (m, q) + T re kn (m, q) (24

)
where T k n is the Fourier coefficients of the convected normal derivative T k . Became different between the singular and regular behavior. It appears to be efficient to treat the two parts in a different fashion. The modal functions (g k n , T re kn ) is the Fourier coefficients of the regular functions (g k , T re k ) are evaluated numerically with standard Gauss-Legendre quadrature, and the singular terms of the static modal Green function G 0 n and its modal convected normal derivative T si kn will be evaluated analytically. Substituting the normal distance in Eq. ( 21), and the using of convected radius Eq. ( 8) of the 2π-periodic even function, the static modal functions (G 0 n , T si kn = d k nq (G 0 )) or the Fourier coefficients of the singular functions (G 0 , T si k = d k n Q (G 0 )) given by Eqs. ( 19) and (21) can be expressed in the form

G 0 n (m, q) = 1 + α 8πcτ * b 1/2 n (α) (25)
T si kn (m, q) = -r nπ/2 t 0 n + n rq r m 2 t 0 n+1 (m, q) + t 0 n-1 (m, q) + ik ε M ∞n G 0 n (m, q) (26)

where the modal function t 0 n is the Fourier coefficient of the function t 0 = 1/4cR * 3 and of the azimuth angle β satisfies the following relation

t 0 n (m, q) = (1 + α) 3 8πcτ * 3 b 3/2 n (α) (27)
The function b s n (α) define the Laplace coefficients of argument α and of order n, the symbol s means a half odd number. These Fourier coefficients can be defined by [START_REF] Izsak | Laplace Coefficients and Their Newcomb Derivatives[END_REF] b

s n (α) = 4 π π 2 0 cos(2nθ) (1 -2α cos(2θ) + α 2 ) s dθ (28)
The argument α of Laplace coefficients in the integrand of the static modal Green function and the modal function t 0 n , can be written as [START_REF] Izsak | Laplace Coefficients and Their Newcomb Derivatives[END_REF]

α = τ 2 (1 + √ 1 -τ 2 ) 2 (29)
The coefficients and their derivatives with respect to α for s = (1/2, 3/2), are all linear functions of the complete elliptic integrals with argument α explained by a particular of the first and second coefficients [START_REF] Brouwer | Methods of Celestial Mechanics[END_REF] b

1/2 0 (α) = 4 π K(α) and b 1/2 1 (α) = b 1/2 -1 (α) = 4 π K(α) -E(α) α (30) b 3/2 0 (α) = 4 π (α 2 -1)K(α) + 2E(α) (1 -α 2 ) 2 and b 3/2 1 (α) = b 3/2 -1 (α) = 4 π (α 2 -1)K(α) + (α 2 + 1)E(α) α(1 -α 2 ) 2 (31)
where K(α) and E(α) are the complete elliptic integrals of the first and second kinds in which the upper limit is π/2 will be respectively denoted by [START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF] K(α) =

π 2 0 1 1 -α 2 cos 2 (θ) dθ and E(α) = π 2 0 1 -α 2 cos 2 (θ)dθ (32)
The difficulty of the analytic evaluation of the static modal Green function and its derivative for largest value of order n is the necessity to develop a recursive relationship between these Fourier coefficients. Using the relations of Laplace coefficients and their derivatives developed in Celestial mechanics calculations [START_REF] Brouwer | Methods of Celestial Mechanics[END_REF], the recursive formulas take the form

b s n+1 (α) = n(1 + α 2 ) (n + 1 -s)α b s n (α) - (n -1 + s) (n + 1 -s)α b s n-1 (α) (33)
The static modal functions G 0 n and t 0 n given by Eqs. ( 25) and ( 27) depending directly of the Laplace coefficients. Then, they are satisfied the recursive formulas Eq. ( 33), which are used the first and the second Laplace coefficients Eqs. ( 30) and (31) to calculate the first terms of these functions. We observe that the static modal convected function T si kn is expressed by an extremely simple formula written in convection term depending of the static modal Green function G 0 n and the modal function term t 0 n evaluated by a simple analytical recursive formula. For the axisymmetric modal case (n=0), the axisymmetric static convected Greens function G 0 0 and its associated convected derivative T si k0 are evaluated by the complete elliptic integrals of the first and second kinds, respectively, K(α) and E(α) of the first and second Laplace coefficients Eqs. ( 30) and (31), can be used to evaluate and simplifies the convection terms [START_REF] Zhang | A coupled FEM/BEM formulation for acoustic radiation in a subsonic non-uniform flow[END_REF][START_REF] Beldi | Some new results for the study of acoustic radiation within a uniform subsonic flow using boundary integral method[END_REF].

This method determined the static modal functions (G 0 n , T si kn ) reduced the complexity of the conventional approximate solution for the static Greens function and its derivatives, the normal derivative and the flow direction derivative and not contain the supplementary terms of the method binomial expansion [START_REF] Juhl | Non-axisymmetric acoustics propagation in and radiation from lined ducts in a subsonic uniform mean flow[END_REF], the simplification of the fast Fourier transform techniques [START_REF] Kuijpers | An improved acoustic Fourier boundary element method formulation using fast Fourier transform integration[END_REF]. When the source point m is taken on the generator with m → q (τ → 1) corresponding of the limit argument α → τ 2 , this has reduced the relative error between the analytical methods used the parameter τ and this analytical method of the Laplace coefficients depending of the argument α, another way between the complete elliptic integrals of the parameter τ and of argument α.

Analytical evaluation of the singular integral terms

Using the convected radiation conditions and the classical procedure to isolate the singularity parts of modal Green function and its convected normal derivative Eqs. ( 23) and (24), we get the new original form of the integral representation in a uniform flow for a point m of the acoustic domain

c ± (m) 2π p n (m) = ± Γ p n (q)T k n (m, q)r q dΓ q ∓ Γ G k n (m, q)d nq p n (q) r q dΓ q (34)
where c ± (m) is the convected angle related to singular integrals derived from the acoustic exterior and anterior problem has been determined by evaluation of the following simple integral, can be rewritten as

c ± (m) = 1 ± 2π Γ T si 00 (m, q)r q dΓ q (35)
The Integral term of convected angle Eq. ( 35) multiply by the revolution angle cylindrical 2π interpreted as the convected angles at the center m. For the point m in an acoustic exterior domain c + (m) = 1 and in the acoustic interior domain c -(m) = 1. When m a regular point on the generator c ± (m) = 1/2. In order to neglect the axisymmetric acoustic pressure contribution at infinity and used the exclusion procedure of the singularity, we find the same form of the axisymmetric BEM developed for Poisson equation in static case [START_REF] Lennon | Boundary integral equation solution to axisymmetric potential flows[END_REF] and Laplace equation in no-flow case [START_REF] Kuijpers | An improved acoustic Fourier boundary element method formulation using fast Fourier transform integration[END_REF][START_REF] Seybert | A special integral equation formulation for acoustic radiation and scattering for axisymmetric bodies and boundary conditions[END_REF] which reduces the computational burden of conventional methods in a subsonic uniform flow and then becomes execution times of no-flow case.

Axisymmetric cylindrical duct

Consider an infinite axisymmetric cylindrical duct of truncated domain Ω i of radius R = 0.1m, of length L = 0.5m, of symmetry axis (oz) and the revolution angle θ, of flow Mach vector M ∞i = M ∞i e z limited by the rigid wall Γ R of normal vector n Γ R and two generators, the left Γ -in z = 0 and the right Γ + en z = L of normal vector respectively n Γ -and n Γ + (Figure 3).

Analytical solution

The axisymmetric acoustic propagation represent by a modal acoustic pressure (n,n') of the azimuthal order n and the radial order n', in the axisymmetric duct checking the axisymmetric wave equation into the harmonic regime Eq. ( 4) is given by [START_REF]Acoustics of ducts and muffler[END_REF] p

± n,n = A n,n J n (k r,n,n r)exp(ik ± z,n,n z) (36)
where A n,n is the axisymmetric normalization amplitude, k ± z,n,n and k r,n,n are the axial and radial wavenumber (± indicates the propagation mode in the increased and decreased z-direction), J n is the Bessel function of the first kind of order n. Since, the term depending of tangential derivative function σ t over the axisymmetric generator duct is zero. Thus, the new condition of radiation or Neumann condition given by the particular derivative of the acoustic pressure Eq. ( 15) where right generator Γ + , resulted on

σ + = ic 2 k ± z,n,n p (37)
We note k ± n,n = c 2 k ± z,n,n the new module wavenumber lower to that of k ± z,n,n . Indeed, when the increasing flow, the axial wavenumber module |k ± z,n,n | is increased but the new axial wavenumber module |k ± n,n | will shrink to a zero value. This explains the efficacy of the new radiation condition for cylindrical duct in a subsonic uniform flow that is tend before the classical radiation condition σ c = ik ± z,n,n p, within the no-Flow case, this condition returns to the classical case and replace the particular normal derivative by the normal derivative (see Figures 4 and5).

Numerical implementation

The numerical implementation of the Axisymmetric Boundary element method Eq. (34) requires the discretization of the interior axisymmetric duct surface Ω i by triangular elements structured and regular with 3 nodes isoparametric of two straight perpendiculars. Therefore, the correspondante generator Γ = Γ -∪ Γ R ∪ Γ + is divided in N nodes and N-1 elements of the convergence criterion in a subsonic uniform flow explain by the relation n * e h * e ≤ λ where n * e is the convected number of finite elements for a radius and axis wavelength λ and h * e is the convected diameter of the finite element, into taking account convected angle in point m = Γ -∪ Γ R and in point m = Γ R ∪ Γ + of the axisymmetric duct c -(m) = 3/4 (Figure 6). The elementary integrals related to the modal Greens function and to the convected normal derivative function of this approximate boundary integral equation in a collocation point m are evaluated analytically by the recursive formula and numerically by Gauss-Legendre quadrature standard method which uses the Acoustics group codes for numerical resolution.

Axisymmetric mode

For the axisymmetric mode of n=0 azimuthal order, we apply the analytical solution (S.Anal) associated to the axisymmetric acoustic pressure without and with flow propagated in the increased z-direction P + n = P + 0,n , BEMA with these new operators, and the FEMA [START_REF] Reddy | An introduction to finite element methods[END_REF][START_REF] Taktak | Numerical Modelling of the Acoustic Pressure Inside an Axisymmetric Lined Flow Duct[END_REF] as a function of the new non-reflection condition Eq. ( 37) where the elementary integrals evaluated analytically because the integrands are regular polynomials except that the second member will be evaluated by Gauss quadrature standard method since it depends on the Bessel function of the first kind. Figures 7 and8 represents the real parts of the axisymmetric acoustic pressure (P + n /A n ) for different values of Mach vector by these three methods. Figures 7 and8 present a real parts of the axisymmetric acoustic pressure (P + n /A n ) by FEMA and BEMA which are in very good concord with analytical solution of relative error given by the L 2 norm. For the axisymmetric finite element method, the relative error is less than 0.5% and for the axisymmetric boundary integral method, the error is less than 0.8%. This error difference between the two methods explained by evaluating the integrals terms of BEMA in the origin point and in the other point of (r=0,z=L). we observed the initial mode (0,0) for the axisymmetric acoustic pressure is uniform along radial and axial wave. The wavelength increases in the case where the increases uniform flow. These explained of the offset wave in the increased z-direction (see Figure 7). In the other mode of axisymmetric acoustic pressure (0,2) from high wave number is non-uniform in a presence and in the absence of flow, along the radial axis and his symmetry axis. Indeed, in the upstream region, the acoustic pressure field is radially amplified to be greater than that obtained in the downstream region. This explains of the Sommerfeld radiation condition Eq. ( 10) and ( 12), and we notice also that the uniform flow is proportional to the wavelength (see Figure 8).

Non-Axisymmetric mode

The second test case concerns a non-axisymmetric mode or a non-axisymmetric excitation represent by the (n = 0, n ) order propagated in the increased z-direction of condition wave ck r,n,n < k i with the cutoff parameter of these mode α r,n,n = k 2 i -ck 2 r,n,n . The figures 9 and 10 shows the modal acoustic pressure obtained by application of the present method such as S.anal and S.BEMA from low wavenumber k i = 40 and a high wave number k i = 80 without flow of Mach vector M ∞ = 0 and with a mean flow of Mach vector M ∞ = 0.4e z . The increase of the subsonic uniform flow from the same mode in the z-direction generates a increase of the wavelength, but for the increase azimuthal order and radial order generates a decrease of the wavelength.

Conclusion

A simple and effective new boundary element method formula is proposed for the axisymmetric acoustic exterior domain and the axisymmetric acoustic interior domain, with or without an arbitrary axisymmetric mean flow, as a function of the two new operators such as particular and convected terms substantially reduce the convective effect of mean flow in the classical axisymmetric boundary integral formulations for original acoustic medium. The main advantage of this method is that the singularity behavior avoids the supplementary evaluation of the modal Green's function and its derivatives like the normal derivative and the flow direction derivative, and requires only the use the simple recursive formula for the static modal Green's function and its convected derivative based on the Laplace coefficients, which reduces the analysis time and the implementation numerical time compared to the conventional formulations and the finite element method.

Figure 1 :

 1 Figure 1: Axisymmetric acoustic problem.

Figure 2 :

 2 Figure 2: Axisymmetric body.

Figure 3 :

 3 Figure 3: Axisymmetric cylindrical duct.

Figure 4 :

 4 Figure 4: |k ± z,n,n (µ n,n )| and |k ± n,n (µ n,n )| for azimuthal order n=0 and the radial order 0 ≤ n ≤ 20 when the wavenumber k i = 5, for different values of Mach number.

Figure 5 :

 5 Figure 5: |k ± z,n,n (µ n,n )| and |k ± n,n (µ n,n )| for azimuthal order n=0 and the radial order 0 ≤ n ≤ 20 when the wavenumber k i = 100, for different values of Mach number.

Figure 6 :

 6 Figure 6: A discretization of the axisymmetric duct surface.

Figure 7 :

 7 Figure 7: Real part of the axisymmetric acoustic pressure (P + 0 /A 0 ) by S.Anal (top), FEMA (middle) and BEMA (bottom) for M ∞ = 0 (left) and M ∞ = 0.4e z (right) when k i = 30.

Figure 8 :

 8 Figure 8: Real part of the axisymmetric acoustic pressure (P + 2 /A 2 ) by S.Anal (top), FEMA (middle) and BEMA (bottom) for M ∞ = 0 (left) and M ∞ = 0.4e z (right) when k i = 90.

Figure 9 :

 9 Figure 9: Real part of the axisymmetric acoustic pressure (P + 1,1 /A 1,1 ) by S.Anal (top) and BEMA (bottom) for M ∞ = 0 (left) and M ∞ = 0.4e z (right) when k i = 40.

Figure 10 :

 10 Figure 10: Real part of the axisymmetric acoustic pressure (P + 2,2 /A 2,2 ) by S.Anal (top) and BEMA (bottom) for M ∞ = 0 (left) and M ∞ = 0.4e z (right) when k i = 80.