Mirlind Bruqi
email: mirlind.bruci@uni-pr.edu

Ramë Likaj
email: rame.likaj@uni-pr.edu

Jorgaq Kaqani

Optimization of Flexible Assembly Systems for Electrical Motors

Keywords: Flexible assembly, optimization, simulation, analysis software, cryptography

The modern assembly systems are fully integrated production systems. They consist of interconnected components of hardware and software. Large flexible assembly systems are being composed of high number of components and managing requires an enormous volume of data. For that reason the controlling and optimization of such systems is too complex. Nowadays that can achieved, using just computers and software tools. In this paper the optimization of product mix for assembly system for electric motors will be described. Also special attention will be given to data security of the information system using cryptography.

INTRODUCTION

Simulation as a modern concept for system analysis, in particular for production systems is driven by rapid technological development of IT. Various experts have developed numerous algorithms, implementing them in various programming languages. Nowadays an engineer's task is to teach the simulation philosophy in order to be set free from very complex mathematical expressions, creating more space for individual creativity and engineering. However for simulation of processes, other knowledge is required depending on the field or the level of research. If we are in a level of a Decision Maker, which requires decisions to be made concerning manufacturing decisions etc, the decision maker should have knowledge of statistics, theory of probability, and a number of other theories linked to buffering and services. Fortunately, modern simulation programs are very advanced, and integrated program packages for statistical preparation of production data and automatic optimization of entire process are included in them. In many cases such simulation process for production system allows very fast and efficient analysis of the stability for the selection of its parameters. Arena is a simulator, which enables fast and efficient modelling of wide natural and artificial systems. It is a language which is designed for analysis and simulation of a large number of production and logistic systems, warehouses and all service systems. Service systems include: restaurants, post offices, banks, hospitals etc. Systems which can be analysed by Arena are listed below:

• Detailed analysis of each type of production systems, transportation devices and industrial robots, • Analysis of complex service systems and customer relationship management systems, • Global analysis of a Supply Chain which include, warehousing, transportation and logistic, • Prediction of system performance based on key parameters such as; cost, quantity, cycle timing, and system efficiency, • Identification of bottlenecks of the Queuing systems and over exploitation of system resources, • Resource planning of staff, equipments and facilities, materials etc. By using Arena and its template the following can be achieved:

• Modelling of relevant processes within imaginary system boundaries, • Simulation of system in order to understand and interpret complex relationship between its elements, identification of strategic points of production segments,

• Visualisation of system operations by graphic animation in order to facilitate both simulation steps, verification and validation, • Analysis of the current state of the system "as it is" or "how it might be", in order to make better planning for the future.

SYSTEM ANALYSIS AND GOAL DEFINING, THE AIM OF SIMULATION

Activities at the flexible system of the electrical machines assembling are very complex, and need to be described in detail. For the purposes of modelling, only main points of them will be subject of analysis.

The system which is considered for analysis from the entry at the first conveyor as a resource point up to the exit from conveyor as an absorbing point for the parts, and is built by following resources:

-Human resources, -Industrial robots, -Pallets with ID chips, -Rolling conveyors. In Figure 1 is shown graphically the path of all entities that circulate in the flexible assembling system. The path of material flux through working stations contains following elements:

•
Input / Output points of the system given by rectangle,

•
Stations given by small circle denoted by R (Resource), • Logic points or system nodes given as parallelepiped (N-Node),

•
Straight lines represent the conveyors, •

Arrows show the direction of the material flux.

Fig. 2. Scheme of material flow through the assembling flexible system

System input represents the source of entities, indicating the input boundary of the system, while output or so-called disappearance point of entities is output boundary of the system. But, the subject of study is analysis of the inner part of the system by considering that all external activities are reduced in two points, input and output. Such approximation brings to the definition the system boundaries.

System allows change the number of resources at working station and also change of material flux depending on overfeed of the system with pallets and to change the frequency of material flux. Therefore, the system input and number of the engaged resources are considered as variables. In the other side, the system is influenced by different demands of the consumers which are made in different timing and in various quantities. These are representing stochastic effects that interfere and put condition to the system itself. The aim of the simulation is to achieve a maximum level of rent ability based on these demands, conditions, schedule and resources.

It is adopted that: C Ri -are the costs for i resources for the simulation period, i=1,..,n, C Mi -are the material costs of specified product i, for i=1 to n (Euro) C Ti -are the costs for waiting time of product-Entity Flow Time i of the system, for i=1 to n (Euro)

Ψ -is the objective function.

By the simulation process we are trying to minimise the objective function. This means that if decrease of the costs is achieved the main goal of optimisation is met, meaning that "with minimum utilisation of resources the maximum production volume is obtained, satisfying technical and organization conditions of the system".

Mathematical model is expressed by:

min( Ψ )= min [ ∑ ∑ ∑ = = = + + n i n i n i C C 1 Ti 1 Mi 1 Ri C ] … (1)
If k-simulations are executed then function (1) has form as given in (2):

min( ∑ = Ψ k j 1 j )= min    ∑ = k j 1 [    + + ∑ ∑ ∑ = = = n i n i n i C C 1 Tij 1 Mij 1 Rij ] C …(2)
So, the simulation that in the best way fulfils the criteria is adopted. This operation is realized with OptQuest package which will be described later.

At OptQuest is chosen the option:

min( ∑ = Ψ k j 1 j ) ≈ = minimize{Entity.Flow.Time} … (3)
In this case, the only goal was the definition of the objective for our experiments in the assembling flexible system. A flexible assembly system is an extremely complex production system. This attribute of the system can be described as the reaction of the market, which forces the manufacturer for decreasing a time to bring up the variety products to market. To accomplish this task, the flexible assembly system must quickly respond to the needs of the market. In this case, the system must be able to reconfigure his structure very quickly. These requirements have a large influence on the organization and technology of a production. Effect of these demands is constantly growing complexity and forces the manufactures to strengthen the flexibility of manufacturing and assembling systems.

To achieve these goals, systems must accomplish the following tasks:

• Increase productivity.

• Increase the quality of products and handle integrated quality management with higher efficiency. • Higher utilization of machines, and stations.

• Quick respond on the external influences and demands. In order to achieve these tasks and to reduce the complexity of the system itself, the complex global objectives can be transformed into a many simple solutions, which are able to meet local objectives. These local objectives are based on hierarchical decomposition of time horizon. In this hierarchical framework there are some interconnected decision levels:

• Flexible assembly system for design (long term).

• Flexible assembly system for planning (medium term).

• Flexible assembly system for scheduling and control (short term). The long term of time horizon means, the period of time on which the decision for design of system is long enough. In that case, in order to reduce the volume of data and the complexity of design process, a lot of activities will be integrated in single assumed activity. The short term of time horizon means, the period of time by which the decision about system and its components will be made. This short term of time horizon give us an answer about the allocation of system resources in a very short period of time. In this paper system is presented the simulation model for electrical motors assembly. This model gives the answer about system behavior. Simulation is made in simulation package ARENA. With OptQuest for Arena is carried out the optimization of product mix.

DESCRIPTION OF FLEXIBLE ASSEMBLY SYSTEM

The core system consists of industrial robots, assembly stations, transport systems and assembly palettes with ID and memory component. The Assembly stations are thought to be automatic transport systems with few integrated assembly sub-systems. Assembly systems are closed cycle. The assembly pallet is the basic carrier for a product during the assembly cycle. The assembly pallets are moved by conveyer belt system. After each assembly step the quality control check is done. If the state of quality check is positive, then the next assembly operation takes place. If the state of quality check is negative, then palette has to go to the repair station and wait for the shop floor operator. Shop floor operator must try to repair the product. If this is not possible he/she has to remove its content from the system, reset the empty pallet and send it in the direction of loading/unloading station. If repair is done, operator must proceed with palette ID and send it to the assembly cycle. New ID corresponds to the current state of product. Palette with repaired product will be rejected and move to the correct assembly station which is defined by given ID. After completing the assembly cycle, palette must be in a loading/unloading station. After that, finished product will be unloaded and palette will be given a new status. This system is shown in fig. 3.

ENTERING SAMPLES FOR FLEXIBLE SYSTEM FOR ASSEMBLING

Another important point of the simulation is to generate data for the working stations known as samples, which is done in a manual and automatic way. Once the system is equipped by a device for identification of pallets and because the system has the possibility of communication with central computer, in this way there is possibility for automatic tracking of data. The data that must be entered in this case are listed in the following:

• Time of arrival of piece at the x station and type, • Time of piece processing at the station, • The time between the fall of assembling system and the duration of the fall, • Schedule of the system, short and long breaks, changes and other organisational restrictions.

BASIC CHARACTERISTICS OF FLEXIBLE ASSEMBLY SYSTEMS

The Simulated flexible complex assembly system deals with following sources of uncertainty and dynamic events:

• Variable availability, arrival rates and has the incoming parts.

• Variable quality of incoming parts -defective parts have to be removed before, either causing errors in system, or being assembled into finished products.

• Mixed batch of parts to be assembled -The assembly station assembling a range of variants of electrical motors types. • Variable availability of resources -resources are tools and machines by which the assembly station must be coordinated directly. • Response to dynamic events in assembly system -For example opening alternative assembly stations to bypass machine breakdowns or allowing more palettes to increase production capacity.

OPTIMISATION USING OPTQUEST FOR ARENA

OptQuest for Arena is designed for manufacturing or business process consultants, analysts and industrial or systems engineers. It is typically deployed as an enterprise business analysis and productivity tool. OptQuest enhances the analysis capabilities of Arena allowing searching for optimal solutions within simulation models. In order to optimize the flexible assembly system in a case study, first the parameters have to be selected which have the largest influence on the product costs. In this case the system is very sensitive to these two parameters:

1. The number of operators in stations, and 2. The number of products on the batch mixes.

Constrains and profits for each Motor type are resented in Table 1: The numbers of resources in different stations may vary, as shown in Table 2.

OPTIMISATION RESULTS

This paper presents simulations model of existing complex FAS. Simulation model consists of 15 different stations. Additionally, quality control of the system is performed after each assembly station. Purpose of the simulation model is to find the best product mix and optimal numbers of replications on OptQuest for Arena, results are acquired and presented in the following tables (table 3 and4). System security is achieved by using cryptographic methods, by using an algorithm which is designed specifically for FAS. Figure 4 represents an interface for system connection.

CONCLUSIONS

Flexible systems for assembling are part of the most sophisticated systems of the time. As such they have the possibility to pass from one assembly program to another in a very short period of time. The flexibility of such systems is its main property. In this way the Passover from one assembly system to another is made in: automatic, semi-automatic way or manually, by replacing the auxiliary devices of robots and machines. Besides the positive side, these systems also have their negative side. Very often different products are mounted to these systems, and from managerial point of view this is heavy workload because of prior planning and preparations. These preparation tasks include: choosing an optimal or suboptimal scenario of material, energy and information flow, including human resource organization.

The selection of such scenarios is made through different optimization algorithms such as: finite, limited and random enumeration. Module varies from one system to another, thus this module must be built for each specific system. The aim of this Paper is to introduce the OptQuest toolbox as new feature of Arena environment for the optimization of complex production systems. The algorithms for the optimization are described in Rockwell Automation handbook OptQuest for Arena. This paper describes on of the easiest ways for optimization of flexible assembly systems using the OptQuest for Arena. The Results of this study are presented in tables 3 and 4.

Fig 1 .

 1 Fig 1. Flexible system of electro-motor's assembling

Fig. 3 .

 3 Fig. 3. Flexible assembly system

Fig. 4 .

 4 Fig. 4. Interface FAS connection The application enables the encryption of various types of files such as "Doc, pdf, doe, dwg etc". The application features are listed in the following: • Symmetric encryption 256bit • Usage of a non-public hash key • Fast encryption/decryption of large files • Intel x86 platform based on Windows Operating System • Non-public algorithm

Table 1 .

 1 Constrains and profits for each Motor type

		Product Type
		A	B	C	D	E
	Constraints	≥ 13	≥ 10	≥ 12	≥ 25	≥ 30
	[% in mix]
	Costs	17	13	15	28	32
	[Euro / Part]

Table 2 .

 2 Variable number of resources per stationThese constraints can be added very easily to OptQuest for Arena and the optimization process can begin.

	Stations		S6	S10	S11	S12	Repair Line
	Number	of	2 to 3	1 to 2	2 to 3	1 to 2	3 to 4
	Resources

Table 3 .

 3 The optimal number of resources on stations

	Stations	S6	S10	S11	S12	Repair Line
	Number of Resources	2	2	3	1	3

Table 4 .

 4 The optimal number of products on batch mixes.

	Product Type on batch mix
	A	B	C	D	E
	15	10	12	33	30


Eric Benoit

Laurent Foulloy

Fuzzy Nominal Scales

Keywords: Measurement theory, Scales, Nominal scales, Fuzzy Scales

INTRODUCTION 1.Nominal scales

Let us present the general scale formalism in order to introduce the aim of the paper. Any scale can be defined with a symbolism [START_REF] Finkelstein | Representation by symbol systems as an extension of the concept of measurement[END_REF]:

C =<X, S, M, R X , R S , F>( 1 )
-w h e r e X refers to a set of object states, and R X is a set of relations on X. -S refers to a set of information entities [START_REF] Mari | The duality of sensing and coding as a fundamental topic of information science[END_REF], and R S is a set of relations on S. -M, called the representation, is a mapping from X to S.

-F is a one-to-one mapping with domain R X and range R S .

In the case of nominal scales [START_REF] Watanabe | A study of nominal scale[END_REF], R S contains only the equality relation on S,a n dR X contains an equivalence relation on X denoted ~. F is then simply defined by .

(

) 2
The equivalence relation ~defines a partition of the set of object states, and M associates each item of the partition to an information entity. The general symbolism of a nominal scale is:

C =<X, S, M,{~},{=},{(~,=)}> (3)
The easiest way for the definition of a nominal scale is to define a partition and a one-to-one mapping from the partition to the set of information entities.

Depending on the measurement aim, the nature of information entities differs:

-An information entity can be a symbol. This is the simplest case used for example with color measurement. -An information entity can be a subset of symbols. For example, a subset of predicates can be used as a representation [START_REF] Watanabe | Theory of classification to reinterpret the nominal scale for our better understanding of measurement information[END_REF].

In this paper we consider the mechanism that represents each object with a fuzzy subset of symbols. This mechanism, called fuzzy description and introduced few years after the fuzzy subset theory in [START_REF] Zadeh | Quantitative fuzzy semantics[END_REF], is now usually applied in fuzzy sensors [6][7].

Fuzzy description

If numerical values are commonly used to represent measurement results, it is now admitted that some applications would better manipulate symbolic values or linguistic terms. For example, a rule-based decision system uses linguistic term in its rules. The advantages of using linguistic terms instead of numeric ones is to reduce the number of symbols.

A measurement result is represented by a fuzzy subset of linguistic terms. The conversion from numerical to linguistic representation is called a linguistic description.

To perform a linguistic measurement, it is necessary to clearly specify the relation between linguistic terms and numbers [START_REF] Benoit | Symbolic sensors[END_REF].

Let X be the universe of discourse associated with the measurement of a particular physical quantity. In order to linguistically characterise any measurement over X,letW be a set of words, representative of the physical phenomenon. Denotes R the relation defined on X×W that formalize the link between items of the universe of discourse and the words of the lexical universe. In order to simplify the management of this link, two mapping are deduced form this relation. The first one is a one-to-one mapping m called the meaning of a linguistic terms [START_REF] Zadeh | Quantitative fuzzy semantics[END_REF]. It associates any word w of W with a subset of X i.e. an item of the set S(X) of subsets of X.

F ~= (,) {} = x i X ∈ () ∀ x j X ∈ () ∀ x i x j ∼ Mx i () ⇔ ,, Mx j () =

Fuzzy nominal scales

Eric BENOIT and Laurent FOULLOY LAMII, ESIA, Université de Savoie, B.P. 806, F-74016 Annecy cedex, France Eric.Benoit@univ-savoie.fr, Laurent.Foulloy@univ-savoie.fr Benoit E., Foulloy L., "Fuzzy Nominal Scales", Proc. of IMEKO TC7 Symposium, Cracow, Poland, 2002, pp. 21-25.

(4)

The other one is a one-to-one mapping d called the description of a measurement. [START_REF] Zadeh | Quantitative fuzzy semantics[END_REF] From the linguistic description point of view, it is interesting to define R such that each object state x is characterized by only one word w: d(x) = {w}. This means that the range of the mapping m is a partition on X i.e. each m(w) is an equivalence class. Thus, the mapping m defines an equivalence relation ~on X. Then the description mechanism can be considered as a nominal scale where: [START_REF] Benoit | A fuzzy colour sensor[END_REF] -The universe of discourse is the set of object states.

-Information entities are singletons of the lexical set W.

-M is the description d.

If the number of information entities is small then the precision of this symbolism is said low. When object states c h a n g ei nas a m ee q u i v a l e n c ec l a s s ,t h ea s s o c i a t e d information entity does not change. Introducing the vagueness into the relation between object states and information entities is a solution for increasing the precision of a symbolism.

Consider now that the relation R is a fuzzy relation on X×W. It is characterized by its membership function denoted µ R (x,w). Then the fuzzy meaning of a word and the fuzzy description of an object state are respectively defined as:

(7) (8)
Where FS(M) is the set of fuzzy subsets of M,andFS(X) is the set of fuzzy subsets of X. These two mappings are simply linked by (9)

As for a crisp description mechanism, it is interesting to define R such that the range of the fuzzy mapping m be a partition on X. Then, a fuzzy equivalence relation ~on X, also called similarity relation, can be deduced. Several fuzzy partitions can be defined. The preferred partitions used for the description with linguistic terms respect the following constraint:

.( 1

THE FUZZY DESCRIPTION MECHANISM IS A SCALE

In this paper, some answers to the following questions are proposed:

-Can it be possible to interpret the fuzzy description of objects into the general scale formalism? -If possible, what kind of scale uses a fuzzy description as representation mechanism?

Considering fuzzy subsets of words as information entities

One solution for the interpretation of the fuzzy description into the scale formalism is to define the representation set as the set of fuzzy subsets of words, and the mapping M as the fuzzy description.

S X () → wW ∈ ∀ mw () xX ∈ xRw {} = dX S W () → xX ∈ ∀ dx () wW ∈ xRw {} = CX w {}wW ∈ () {} d {} = {} ~= , () {} ,, , , , ¢² = mX F S W () → wW ∈ ∀ xX ∈ ∀ ,µ mw () x () µ R xw , () = dW F S X () → wW ∈ ∀ xX ∈ ∀ ,µ dx () w () µ R xw , () = wW ∈ ∀ xX ∈ ∀ ,µ dx () w () µ mw () x () = xX ∈ () ∀µ mw () x () wW ∈ ¦ , 1 = µ ~xy , () min µ mw () x ()µ mw () y () , () () wW ∈ ¦ = {0.8/

Denote MFS(W) the subset of FS(W) such that each item of MFS(W) verifies (14). It's interesting to see that MFS(W)

is the union of ranges of all possible fuzzy descriptions on the lexical set W.

(15)

In this solution the symbolisation is precise. This property come from the size of the set of information entities that can be very large. This is also contradictory with the first interest of fuzzy description i.e. reducing the number of information entities in order to simplify the manipulation of these entities.

Considering fuzzy description as the representation mapping

In an other solution, the information entities are symbols, and M is a fuzzy relation from X to the lexical set W.

(

-w h e r e X refers to the same set of object states, -where information entities are words of the lexical set W, -w h e r e R the fuzzy relation from X to W that defines the fuzzy meaning and the fuzzy description ( 7) and ( 8).

.

Fig. 2. Scale interpretation

In this case, the symbolisation stays precise, because the equivalence relation on object states is the same than in preceding case. The number of representational symbols is small. This is useful when we need to handle the knowledge associated with object states.

It can be noted that using a crisp partition instead of a fuzzy one induces a usual nominal scale. Introducing a fuzzy description of object states instead of a crisp one induces a proximity relation between items of the fuzzy partition and then induces a proximity relation between information entities that represent these objects states. As symbolisms are used in order to give tools for manipulating information entities instead of their associated object states, a symbolism do not have to introduce relations that not exist between object states.

FUZZY NOMINAL SCALES

Using fuzzy description can be an approach for taking into account existing proximity relation between object states. And fuzzy nominal scales can be used when nominal scales are too poor and when an ordinal scale cannot be defined.

COLOR MEASUREMENT EXAMPLE

The color measurement is a typical measurement where its seems to be possible to use a scale more informative than a simple nominal scale. Indeed, we can imagine that a proximity relation can exist into a color space. Red color is close to orange color but not close to green color. Unfortunately no existing crisp scale presently includes this proximity relation.

A color measurement is performed with three photometric transducers that recreate the effects of red, green and blue human cones. The set of color measurements is part of the so-called RGB cube. Then a non linear mapping is applied

µ A w () () wW ∈ ¦ 1 = AB , () MFS W () 2 ∈ ∀ µ ≡ AB , () min µ A w () µ B w () , () wW ∈ ¦ = CX W R {} = {} ~= , () {} ,,, , , ¢² = wz , () W 2 ∈ ∀ w=z w {} ≡z {} ⇔ x1 x2 x3 X M w 1 w 2 w 3 0.8 0.2 0.4 0.6 1 W wz , () ∀ W 2 ∈ µ = wz , ()sup Ad X () ∈ min µ A w () µ A z () , () () = w ∀ W ∈ sup Ad X () ∈ µ A w () () 1 =
to the RGB cube. This mapping has been chosen in order to obtain a colour information which do not depend on the luminosity. This information, called the chrominance, is included into the chrominance plane. In this example, we choose the C1C2 chrominance plane [START_REF] Benoit | A fuzzy colour sensor[END_REF]. Seven typical colors are localized on this plane (Fig. 3.).

Fig. 3. C1C2 chrominance plane used for color measurements

A fuzzy partition based on an interpolation of human knowledge can be defined on this chrominance plane (see Fig. 3.). The human knowledge is defined by a set of measurement states, With each state associated to a word. The interpolation knowledge is performed with the triangulation method. Other methods like fuzzy C-means, or K-Nearest neighbour that can also be used are not presented in this paper. This fuzzy partition is used for the definition of the fuzzy nominal scale. The proximity relation on the lexical set can then be define (see table 1).

The equivalence/proximity relation = on color words shows that red colors are near yellow colors. In can be noted that if the lexical set is larger, including for example the orange color, then the relation = can differs, and red colors can be near orange colors but not near yellow colors. Actually, the definition of the proximity relation on the lexical set of words depends on the definition of the fuzzy partition. It depends on the fuzzy meaning of each word.

CONCLUSION

In this paper, the definition of a fuzzy nominal scale had been presented, and the interest and the danger of using such scale where exposed. Using fuzzy nominal scale for a measurement induces a proximity relation on the representational set that is more informative than the equality. This is an improvement only if this proximity relation represents a real proximity relation between object states. If not, the danger is to introduce a relation that does not have any meaning. As several family of fuzzy equivalence relations can be used, it will be useful to study the properties and the useability of associated fuzzy nominal scales.

  0 ) This constraint reveals the need of exclusion between linguistic terms for the description of states. Based on this partition a similarity relation can be extracted: (11) This relation is reflexive, symmetric and verifies a weak version of transitivity based on the Lukasiewicz T-norm. With this fuzzy equivalence relation it seems possible to interpret the fuzzy description as a nominal scale.

Fig. 1 .

 1 Fig. 1. Scale interpretationThe following symbolism can then be defined:

  mW

  -the relation ~is the fuzzy equivalence relation on X definedin(11), -a n d = is an equivalence relation on items of W.T h i s relation can be simply defined such that:

Fig. 4 .

 4 Fig. 4. Human knowledge and triangulation of the chrominance space.

Fig. 5 .

 5 Fig. 5. Fuzzy meaning of the word Red.

  To the question «what kind of scale uses a fuzzy description as representation mechanism?», the last symbolism presented seems to let conclude that a fuzzy description is a nominal scale. But we consider that this conclusion is dangerous. Indeed, in the case of a nominal scale, the only relation that can be used on the lexical set is an equivalence relation. But the introduction of fuzzyness into the measurement mechanism induces a new proximity relation on the lexical set. This relation comes from relation between the meanings of words i.e. the items of the fuzzy partition. This relation is an equivalence relation that verify a weak transitivity based on Lukasiewicz T-norm. This relation does not verify a strict definition of the transitivity i.e. the sup-min transitivity. It is then considered as a proximity relation. This relation is not the crisp equality relation, and can link different words. Actually it acts as a proximity relation between words.

	As each word w belongs to an item of d(X) it is possible to
	define an equivalence relation based on the general
	combination-conjunction:
	(18)
	such that
	(19)


Mônica Alessandro Da Silva

Franchi Carniello

José Da

Alessandro Luiz Da Silva

Mônica Franchi Carniello
email: monicafcarniello@gmail.com

José Luís Gomes Da Silva
email: gomesdasilvaster@gmail.com

D A Silva

M F | Carniello

J L G Silva

Holistic Vision of Sustainability in the Production Chain in Oil Exploration Pre-Salt Layer

Keywords: Productive Chain, Regional Development, Pre-salt, Sustainability Geographic

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

Introduction

The current global scenario, under sustainability aspects, reflects the thinking reflected in the work of the Club of Rome in 1972, which indicated the reality of the scarcity of natural resources, considering the high rate of world population. With the emergence of "environmental crisis" a contemporary line of thought gained political space. Thus, current discussions about the future, considering, especially the development has been promoting joint socio-political in favor of sustainable growth, with relevance to the thoughts of [START_REF] Sachs | Estratégias de transição para o século XXI -desenvolvimento e meio ambiente[END_REF], which emphasizes the perspectives, social, economic, environmental, space and cultural sustainability, particularly considering the exploration of the Brazilian pre-salt, particularly considering the exploration of the Brazilian pre-salt, observed by outlays on research and technological development, since there will be high demand for technical training in logistics and risks of exploration.

Purpose

The objective of the research is in the discussion of the current global scenario, under sustainability aspects, considering a holistic vision of sustainability in the production chain in oil exploration pre-salt layer which in a way reflects the thinking of the work referenced in the Club of Rome in 1972 which indicated the reality of scarcity of natural resources, considering the high rate of world population. With the emergence of "environmental crisis", a contemporary line of thought gained political space. Thus, the current discussions about the future, considering in particular the development have promoted socio-joint policies for sustainable growth, with relevance to the thoughts of [START_REF] Sachs | Estratégias de transição para o século XXI -desenvolvimento e meio ambiente[END_REF], which emphasizes the perspectives, social, economic, spatial, environmental and cultural sustainability.

Methodology

The methodology of the article provides a qualitative approach, based on literature search, design specification, aiming for a holistic supply chain to discuss the pre-salt layer, from the perspective of sustainability, with relevance to the logistics flow. The research was the period between 2005 and 2010, considering scientific publications that give context to the pre-salt exploration, and the relevance of complexity from their extraction and mainly for transport, given that the sustainability aspects refer to the deep thoughts, considering the of environmental events in recent years. The research was defined as literature, with a descriptive design with a qualitative approach. In addition to scientific publications, to collect data, we also used the database of the Center for Research and Development (CENPES) considering, specially Petrobras's information, relevant to planning and exploration of technological advances in the layer of pre-salt.

Results and Discussion

As a result, it was found that since its discovery in 2006, the exploration of the pre-salt layer is highly complex, requiring high expenditure on research and technological development, especially in innovation, reaching its fullest designed for 12 years, considering the problems under the aspects of environmental, cultural, economic, and social space. Considering the extraction of oil in the country, with relevance to the pre-salt layer since its discovery in 2006, provides a holistic view of sustainability and understanding of its complexity, in particular its operation, but is concentrated phase of research and development. The development maximized and designed for 12 years since its discovery in 2006, certainly, will require further large investments, mainly in Research, Development and Innovation (RD & I), given the problems in aspects of environmental, cultural, economic, spatial and social considering the complexity and magnitude of the presalt exploration, surrounded by risks, among which we highlight the sustainability space, reflected mainly by logistical problems, considering the geography of the pre-salt, and environmental impacts involved in this process extraction, these impacts still purely theoretical knowledge, as represented at Table 1, by stages since company foundation.

Years Since Petrobras Foundation (1954)

Goal of Production Thus a holistic view represented in Fig. 1, under a projection of an economic standpoint, the pre-salt will be explored when likely to be high oil prices, yet in a not too distant future, low oil prices, due to the heating global or the development of cleaner sources of energy, can impair the operation of the pre-salt.

Considering the array of sustainability, Lima and Fernandes (2009) a holistic view, make an approach on the oil sector, where the exploration and development stages express the relevance of the processes detailed below:

• Exploration and delineation of mineral deposits; Drilling, completion and production activities are concentrated in drilling the well and prepare it to receive the permanent production equipment for oil and gas;

• Refining and transportation back to the transfer of crude oil and gas production units for the derivatives, while the processing and distribution of natural gas cover asset of operations for the transport;

• Distribution and use of natural gas and other hydrocarbon fluids.

Geography of the Brazilian Pre-Salt Layer

Geographically the pre-salt layer, covers part of the north coast of São Paulo, with a huge reservoir of oil and natural gas, also present in the Santos, Campos and Espirito Santo, and presents a macro view with coverage in the coastal region between the states of Santa Catarina and the Espirito Santo. Above all, the challenges of the pre-salt exploration beyond the technological and environmental aspects, with the geographic relevance to sustainability, given the logistical problems of transport, considering the existence of a huge volume of natural gas, inherent in the process. Variable important, observed particularly in the research of [START_REF] Lima | Os desafios, os impactos e a gestão da exploração do pré-sal[END_REF] which notes that there is a huge gap in about 350 km the continent almost making it impossible to flow of gas through pipelines across the ocean, and thus it is necessary to burn it at the base of production. So the question is:

• Considering sustainability, which the environmental impacts of burning of these gases, as it focuses on issues related to the high release of gases?

• The same concept of sustainability with geographical approach, how will the distribution logistics of this production so far inland from the coast?

Thus, this paper discusses the holistic view of sustainability, considering the oil exploration of the pre-salt layer, in particular from the perspective of geography.

Geography of Global Pre-Salt Layer

According to a survey of Lima (2008, p. 4), it is estimated that "the pre-salt layer has an area of 112,000 km ². Of this total, 41.000 km ² has been tendered and awarded. This layer is about 800 km long and in some areas, 200 km wide, "considering the reserves in Brazil, located at 350 kilometers of coastline. And also, the pre-salt layer is located between five and seven thousand feet below sea level, it can reach up to two thousand feet thick,

Challenges of Sustainable Development of the Environmental and Cultural Vision of Oil Supply Chain

The concept of sustainable development is concatenated from the emergence of the environmental crisis in mid-1980. Currently, "[...] is closely related to action strategies to implement a set of measures and therefore is not just a theoretical concept, but only instrumental". (BARRETO, 2001 cited in Fernandes et al. 2007, p. 15). [START_REF] Sachs | Estratégias de transição para o século XXI -desenvolvimento e meio ambiente[END_REF] describes the cultural sustainability, assumes the role of relevance in the process of regional sustainable development, considering the roots of indigenous models of modernization, especially the regulatory processes of continuous improvement. Considering the approach of Jacobi (1999), "the issue of sustainability is, in this end of century, a central role in the reflection on the dimensions of development and alternatives that are configured" in a way, "[...] the social and environmental context that characterizes contemporary societies reveals that the impact of humans on the environment "[...] is becoming increasingly complex, both in quantitative and qualitative terms." (JACOBI, 1999, p. 175, emphasis added).

In the same context, recent research by Rodrigues (2010, p. 26) states that no doubt "one of the challenges faced by the Brazilian Petroleum Company (Petrobras), the exploration of oil in the pre-salt is the great amount of associated natural gas." [START_REF] Rodrigues | Pré-sal: desafios tecnológicos[END_REF], this same approach, reveals that the gas transported by pipeline, by the ocean, considering the enormous distance, it would be economically unviable, and also contextualizes, in his research that an alternative would be burning even in the very basis of production, where, however, would encounter another aggravating factor, the production of pollutant gases on a large scale, versus the environmental dimension as defined: "The environmental dimension is focused on a set of actions that aim to prevent or mitigate harm arising from the development process." The solution is not in the suppression of development actions, but the adoption of forms of management that are intended to use not predatory the planet's resources. (CARVALHO, 2005, p. 49). Considering the magnitude of the complex technological, and as their own definitions of the Brazilian Petroleum (Petrobras), initially for oil exploration in the pre-salt, several alternatives were suggested techniques in relation to excess gas inherent in the process, among which stands out:

• The technologies of the transformation of natural gas for transportation in lieu of using pipelines, which would involve high costs, economically unfeasible;

• Another important variable defined as a challenge to be considered in the process of extracting the pre-salt, is the cultural sustainability with a focus on training of human resources that are necessary for the production. At the same focus, [START_REF] Fraga | O pré-sal e seus desafios[END_REF], estimated that 110 000 professionals will be needed to attend this demand.

But how to meet this demand for human resources?

The technological absorption capacity can be defined by the observations of [START_REF] Cohen | Absorptive capacity: a new perspective on learning and innovation[END_REF], considering that this concept can be understood since it has the ability to recognize the value of new knowledge, which is extremely necessary that assimilation. This learning process inherent in technological absorptive capacity is of fundamental importance, but there is a learning curve -learning curve, that is, time is a valuable expenditure necessary in the absorption of technological knowledge is considered as another variable in process operation of the pre-salt. In this approach, although innovative capacity of enterprises depends directly from their capital borderline intellectual, and how they are structured, considering a range of variables such as geographic location, its investments in RD & I, this same perspective, Sveiby (1998, p. 24), emphasizes that the so-called "knowledge organizations" are those characterized by having few tangible assets and their intangible assets much more valuable than their tangible assets. And According to Stewart (1998, cited by ANTUNES; MARTINS, 2007, p. 6), "[...] tangible assets contribute much less to the value of your product or service end of that intangible assets." Whereas the pre-salt extraction requires much more than prior knowledge, and requires investment in large scale research, particularly De Negri (2006, p.

108) in their study identified some factors to measure the ability absorption of the companies, reducing the learning curve. This requires also was observed by Da Silva; Oliveira; Silva (2011, p. 4), by Factors determining the absorption capacity of technology companies. Considering the operational know-how, it is clear that Petrobras operates 23% of global production in deep water, being independent of this framework, the exploration of the pre-salt, goes further, and should be considered as new business to a new market, given that alternative fuels increasingly gaining space, not to mention the exploitation of "black gold", the shortage will lead to the same.

Another large and important variable to be considered, without doubt are the environmental challenges and their possible impacts, given recent events, where you can highlight the accident by British Petroleum -BP in mid-2010. The accident with the Deep-water, the platform operated by British Petroleum (BP) in the Gulf of Mexico undoubtedly produced much more than the leak of five thousand barrels of oil per day. In addition to the environmental consequences, with a spot of more than 50 km in length, the accident did the British company losing $ 25 billion in market value and raised some questions, among which are:

• Considering the magnitude of environmental accident, on an international scale, to what extent the operation of drilling in high seas has been affected?

• While environmentalists around the world to question this type of operation offshore, pointed to by experts as one of the safest in the production chain of oil black, accidents such as the Gulf of Mexico in 2010, are increasingly rare?

Considering the problems of distribution logistics, with special emphasis on transport, as found in the survey of [START_REF] Rodrigues | Pré-sal: desafios tecnológicos[END_REF], it is important to note that other studies reveal that large oil spills at sea occurred more frequently and were more harmful the environment during transport of the product and not on his farm, the scene observed that the representation in Table 2.

Fig. 1 .

 1 Fig. 1. Macro view of the productive chain of the pre-salt (Source: Authors, adapted from Lima and Fernandes (2009, p. 323).

Table 1 .

 1 Macro view of oil exploration in the country (Source: Authors adapted from Petrobras -Centre for Research and Development -CENPES 2010).

1.000.000 barrels of oil per day

			Chronology
	45	45	1954
	27	27	1974
	22	22	1984
	12	12	2006

Table 1 .

 1 Ranking of the 10 largest oil spill accidents at sea (Source: The International Tanker Owners Pollution Federation -ITOPF 2010).

Date

associated with extraction and production of oil in the Santos Basin. One proposes the installation of filters made of activated charcoal in the bottom of the sea, for the absorption of oil in case of leaks in vessels or platforms. The other aims to implement a regional referral center for study control and monitoring of aquatic and terrestrial environments, in order to protect the biodiversity of the regions explored.

The Petrobras itself (2011) states that the challenges are to find ways to provide efficient and especially in finding logistics solutions for this energy source more and more distant. However the flow of production is still a topic in depth discussion, given the risks and costs involved, but is already approaching the reinjection of gas into the soil. Gouveia (2010, p. 33) also reveals that, "[...] the logistics organization is another bottleneck: How to support the transport of people and cargo and supply of diesel for operation of the rigs and production platforms? The distance of the pre-salt accumulations from the coast is about 300 km, which is the maximum range of flight autonomy of most helicopters. Petrobras's own research shows that the solution may be the construction of intermediate bases, installed between the point of production and shipment of resources, providing better techno-economic feasibility for the process of extracting the pre-salt. However the search for optimal solution for the scenario of the pre-salt, Petrobras has signed partnerships with institutes and universities stand out among the fifteen Brazilian universities:

• Universidade Estadual Paulista (UNESP),

• São Paulo University (USP),

• State University of Campinas (UNICAMP),

• Universidade Federal de São Carlos (UFSCar),

• Federal University of Rio de Janeiro (UFRJ),

• Universidade Federal da Bahia (UFBA),

• Institute of Technological Research of São Paulo (IPT), among others.

Practical Implications

This study aims to contribute about thinking in logistics and transportation of petroleum considering the extraction of the same in the pre-salt layer, where large distances must be considered in order to avoid impacts, particularly environmental.

Value of the Paper

The issue of sustainability, in particular the geographical, directly affects the operation of the pre-salt, considering especially the distances involved, the logistics of transportation and distribution of products, and especially the disposal or reuse of excess natural gas, which is inherent extraction, which somehow has encouraged large-scale surveys in order to obtain effective responses, and provides technical-economic in exploring the viability of the pre-salt layer.

An aggravating factor seen in the research, especially considering the formation of intellectual capital, will be the peak of production that require a highly qualified man power, considering that there will be a learning curve, impacting the absorption curve of technology. How is the timing between research and development, innovation and application themselves of the resources involved? Considering the presence of foreign corporations in the process of exploitation of the presalt layer, as is being conducted to integrate technology, considering the learning curve? Given the magnitude of the project, current events and oil spill at sea. Considering all the variables considered in this and other research, the question is: Are we truly prepared? Thus, the main result was observed in research that somehow the oil and natural gas should and must have a holistic vision for the future of energy, considering the scarcity of resources and innovation in the search for sustainable alternatives.

In this approach, the coming decades are likely to be the most challenging energy history of the planet, especially considering the recent scientific studies that point to alternative fuels, in a more sustainable world in fact.


Hang Zhou
email: hang.zhou@recherche.enac.fr

Alexandre Chabory
email: chabory@recherche.enac.fr

Rémi Douvenot
email: douvenot@recherche.enac.fr.

A 3D Split-Step Fourier Algorithm Based on a Discrete Spectral Representation of the Propagation Equation

Keywords: Atmospheric propagation, discrete Fourier transforms, electromagnetic propagation, modeling, operators, parabolic wave equation, refraction

A propagation method self-consistent in the discrete domain for the simulation of large distance scenarios in a 3D inhomogeneous atmosphere is proposed. First, a spectral representation of the field in a discrete 3D domain over an infinite perfect ground is introduced in cylindrical coordinates. A discrete spectral propagator is then derived. As in 2D methods, the inhomogeneous atmosphere is considered in the spatial domain using phase screens. An explicit numerical scheme is introduced. For directive sources, a sectoral propagation is proposed to reduce the computation time by limiting the computation domain. The method is tested by simulation scenarios involving refractive effects in both vertical and azimuthal directions. The 3D refractive effects are shown to be properly simulated.

I. INTRODUCTION

T HE parabolic equation method is an approximation of the wave equation valid along a paraxial direction neglecting backward propagation. This equation has been used extensively in radio wave propagation [START_REF] Levy | Parabolic Equation Methods for Electromagnetic Wave Propagation[END_REF], e.g., for radar coverage prediction. One of the most efficient method for applying the parabolic equation is based on split-step Fourier methods as pioneered by Hardin and Tappert [START_REF] Hardin | Applications of the split-step Fourier method to the numerical solution of nonlinear and variable coefficient wave equations[END_REF].

These algorithms can take into consideration the antenna pattern, the ground composition, the relief, and the atmosphere. Nevertheless, the backscattering is neglected. In order to evaluate the field propagation iteratively at increasing distances, the computation is realized going back and forth from a spatial to a spectral representation of the wave. Indeed, derivatives become multiplications in the spectral domain which simplifies and alleviates the computations.

The electrical characteristics of the atmosphere are taken into account in the spatial domain using phase screens [START_REF] Levy | Parabolic Equation Methods for Electromagnetic Wave Propagation[END_REF]. The ground is taken into account by choosing a spectral transform consistent with the boundary condition. If the ground is a perfect electric conductor, the spectral transform corresponds to a sine or cosine Fourier transform depending on the considered field. For an impedance ground, the discrete mixed Fourier transform is usually used [START_REF] Dockery | An improved impedance-boundary algorithm for Fourier split-step solutions of the parabolic wave equation[END_REF]. Note also that the ground relief can be accounted by means of a staircase approximation [START_REF] Levy | Parabolic Equation Methods for Electromagnetic Wave Propagation[END_REF]. These methods have initially been introduced for 2D cases, with either an axial symmetry in Cartesian coordinates or a rotational symmetry about the vertical axis in cylindrical coordinates. To model a 3D scene with the latter, a N × 2D model approximation scanning all azimuth angles is classically used [START_REF] Hardin | Applications of the split-step Fourier method to the numerical solution of nonlinear and variable coefficient wave equations[END_REF]. However, N × 2D models neglect lateral effects. These latter can be significant, especially in the case of an irregular relief and/or non-constant refractive index.

In the context of underwater propagation of acoustic waves, 3D finite-difference or split-step Fourier methods have been used [START_REF] Siegmann | A wide-angle threedimensional parabolic wave equation[END_REF] [START_REF] Sturm | On the use of higher-order azimuthal schemes in 3-D PE modeling[END_REF] [6] [START_REF] Fawcett | Modeling three-dimensional propagtion in an oceanic wedge using parabolic equation methods[END_REF]. The 3D models demonstrate a better accuracy than N×2D models where azimuthal couplings are present. The acoustic pressure and electromagnetic fields both satisfy the 3D Helmholtz equation, so the 3D propagation can be developed similarly for acoustic and electromagnetic propagation. Later on, for electromagnetic propagation, an extension of split-step methods based on the parabolic equation to 3D has been developed in Cartesian coordinates [START_REF] Zaporozhets | Bistatic RCS calculations with the vector parabolic equation method[END_REF] [START_REF] Janaswamy | Path loss predictions in the presence of buildings on flat terrain: A 3-D vector parabolic equation approach[END_REF].

The previous 3D methods are based on the parabolic equation, in which a paraxial approximation is applied to split the vertical and azimuthal derivative components. In this paper, this approximation is here bypassed by using an exact spectral representation of the field in a homogeneous atmosphere. This representation is obtained through a diagonalisation of the vertical operator as proposed by Chabory et al. [START_REF] Chabory | An exact spectral representation of the wave equation for propagation over a terrain[END_REF]. Janaswamy has also derived a similar representation [START_REF] Janaswamy | Radio wave propagation over a nonconstant immittance plane[END_REF]. Besides, this operator has been extended to 3D in cylindrical coordinates in [START_REF] Chabory | An exact vectorial spectral representation of the wave equation for propagation over a terrain in 3D[END_REF]. These methods are more accurate since no paraxial approximation is assumed. Indeed, they are derived from the Helmholtz equation, not from its parabolic approximated form.

Since the aforementioned methods are used for numerical simulations, the computed fields or potentials are expressed in the discrete domain. Therefore, a discrete formulation of the propagation equation is derived to achieve self-consistency according to the discrete electromagnetic theory [START_REF] Chew | Electromagnetic theory on a lattrice[END_REF]. The discrete mixed Fourier transform (DMFT) [START_REF] Dockery | An improved impedance-boundary algorithm for Fourier split-step solutions of the parabolic wave equation[END_REF] has been a first effort towards this discrete formulation. Indeed, the boundary condition is given in the discrete domain by expressing the partial derivative along z with a finite-difference approximation. However, the propagator is derived from the continuous spectral equation. This is an inconsistency in the DMFT theory.

In this paper, we fill the gap by deriving a 3D propagator from the discrete equations directly. Consequently, it is selfconsistent in the framework of the discrete electromagnetic theory. The approach is also valid in 2D. The ground and at- mosphere effects are considered in the spatial domain as in 2D methods. For the sake of clarity, the theoretical development is given over a perfect metallic ground. Considering a slowly varying inhomogeneous atmosphere, the 3D split-step Fourier method is applied. Its explicit numerical scheme is introduced.

In Section 2, a homogeneous atmosphere is considered to express the problem in cylindrical coordinates by means of potentials. The propagation is expressed by means of the discrete spectral representation. The discretized form of the 3D formulation is presented. In Section 3, the explicit numerical scheme in an inhomogeneous atmosphere is obtained from phase screens and a 3D split-step algorithm. In Section 4, a sectoral propagation method is introduced. In Section 5, the sectoral propagation method and the 3D formulation are validated by comparison with an analytic expression. Finally, numerical tests are performed in the presence of inhomogeneous atmospheres to show the accuracy of the method.

II. EXACT DISCRETE SPECTRAL REPRESENTATION IN A HOMOGENEOUS ATMOSPHERE

A. Configuration and continuous formulation

We use the cylindrical coordinate system (r, θ, z) with unit vectors (r, θ, ẑ). We assume that the sources are located in the cylinder r ≤ r 0 and that the fields are known at r = r 0 . A perfect metallic and infinite planar ground is located at z = 0. The propagation is computed in the region r > r 0 , z 0.

In a preliminary step, a homogeneous refractive index, i.e. n constant, is considered and the vertical domain is limited to z = z max where a PEC condition is applied. In this case, the fields can be decomposed in one transverse electric (TE) and one transverse magnetic (TM) components with respect to z by means of Hertz potentials oriented along ẑ. For the electric field, omitting the time-dependence e jωt yields

E = k 2 0 n 2 Π e + ∇∇ • Π e -k 0 ζ 0 n∇ × Π h , (1)
where Π e and Π h are the electric and magnetic vectorial potentials, k 0 is the wavenumber in free-space, and ζ 0 is the free-space impedance.

Finally, replacing Π e by Ψ √ r ẑ gives the TM electric field

E = ∂ ∂r r -1 2 ∂ ∂z Ψ r + r -3 2 ∂ 2 Ψ ∂z∂θ θ + r -1 2 ∂ 2 Ψ ∂z 2 + k 2 0 n 2 Ψ ẑ, (2)
where Ψ is the solution of

∂ 2 Ψ ∂r 2 + 1 r 2 ∂ 2 Ψ ∂θ 2 + ∂ 2 Ψ ∂z 2 + k 2 0 n 2 + 1 4r 2 Ψ = 0. (3)
Furthermore, due to the perfectly conducting walls at z = 0 and z = z max , Ψ fulfills the boundary conditions

Ψ| z=0 = 0, Ψ| z=zmax = 0. ( 4
)
Note that the TE case can be formulated in a similar way from a z-directed magnetic potential Π h = Ψ √ r ẑ. The main differences are Neumann boundary conditions and a different expression for the electric field.

In the proposed method, the potentials are propagated. The electric field is obtained afterwards from the potentials.

B. Discrete formulation

For obvious numerical reasons, the computation domain is discretized and of finite size. Thus, the vertical domain is limited to z ∈]0, z max [ and the following uniform grid is used

z = p z ∆z for p z = {1, ..., N z -1}, θ = p θ ∆θ for p θ = {0, ..., N θ -1}, (5)
with ∆z = z max /N z and ∆θ = 2π/N θ . The grid is shown in Fig. 1. To obtain the discretized counterpart of (3), the second-order spatial derivatives with respect to θ and z are discretized on the interior points of the grid using the central-difference approximation. This yields

∂ 2 Ψ ∂r 2 + 1 r 2 d 2 θ Ψ + d 2 z Ψ + k 2 0 n 2 + 1 4r 2 Ψ = 0, (6)
with

d 2 z Ψ pz,p θ = 1 ∆z 2 (Ψ pz+1,p θ -2Ψ pz,p θ + Ψ pz-1,p θ ), d 2 θ Ψ pz,p θ = 1 ∆θ 2 (Ψ pz,p θ +1 -2Ψ pz,p θ + Ψ pz,p θ -1 ), (7)
where Ψ pz,p θ = Ψ(r, p θ ∆θ, p z ∆z).

C. Discrete spectral representation

To render our spectral representation numerically selfconsistent, we directly work from [START_REF] Sturm | Numerical study of broadband sound pulse propagation in three-dimensional oceanic waveguides[END_REF]. Since Dirichlet and periodic boundary conditions are imposed along z and θ, the discrete spectral representation amounts to a discrete sine transform along ẑ and a discrete Fourier transform along θ. Hence, the spectral transform, denoted as Ψ = U Ψ, is explicitly given by

Ψq θ ,qz (r) = 1 √ N θ 1 2(N z + 1) N θ -1 p θ =0 Nz-1 pz=1 Ψ p θ ,pz (r) e -j 2πp θ q θ N θ sin π p z q z N z . (8)
Using fast trigonometric transforms, the computation of either U Ψ or U -1 Ψ can be efficiently performed, typically in N θ N z log 2 N θ N z operations. Applying U to (6) yields

r 2 ∂ 2 Ψ ∂r 2 + (r 2 k 2 r -κ 2 ) Ψ + 1 4 Ψ = 0, (9)
where

κ = 2 ∆θ sin πq θ N θ for q θ = {0, ..., N θ -1}, (10)
and

k 2 r = k 2 0 n 2 -k 2 z with k z = 2 ∆z sin πq z 2N z for q z = {1, ..., N z -1}. (11)
Note that an impedance boundary condition could be considered by replacing the DST by the DMFT [START_REF] Feit | Light propagation in graded-index optical fibers[END_REF].

In a homogeneous atmosphere, the field can be propagated from r 0 to any distance r > r 0 in the spectral domain. Indeed, the analytical solution of ( 9) is [START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF]

Ψ(r) = H (2) κ (k r r) H (2) κ (k r r 0 ) r r 0 Ψ(r 0 ), (12)
where H

(2) κ denotes the Hankel function of the second kind and order κ.

This equation is similar to the propagation of cylindrical harmonics except for the expressions of k r and κ because of the discrete formulation. When q θ N θ , i.e. for modes slowly varying with θ, the order of the Hankel function becomes an integer (κ ≈ q θ ) as in the continuous case [START_REF] Harrington | Time-Harmonic Electromagnetic Fields[END_REF].

For κ is fixed and |k r r| → ∞, the Hankel function can be simplified to reduce the computation load. Doing so, the spectral propagation from r 0 to r becomes

Ψ(r) ≈ exp(-jk r (r -r 0 )) Ψ(r 0 ). (13)

III. EXPLICIT NUMERICAL SCHEME FOR A SLOWLY VARYING REFRACTIVE INDEX

A. Phase screens

We aim at computing the effects of refraction in a 3D configuration in the presence of a planar and perfectly conducting ground. The atmosphere is characterized by a slowly varying refractive index n. To treat an inhomogeneous atmosphere, a method of phase screens similar to [START_REF] Dockery | An improved impedance-boundary algorithm for Fourier split-step solutions of the parabolic wave equation[END_REF] is proposed.

From equation (3), we have

1 k 2 0 ∂ 2 ∂r 2 + (A + B + 1) Ψ = 0, (14)
where

A = 1 k 2 0 1 r 2 ∂ 2 ∂θ 2 + ∂ 2 ∂z 2 + 1 4r 2 , B = n 2 -1. (15)
The equation ( 14) is decomposed as a forward and a backward propagation terms. The backward term is neglected. Therefore, Ψ fulfils

1 k 0 ∂ ∂r -j √ A + B + 1 Ψ = 0. ( 16
)
This formulation is exact for the forward propagation in a homogeneous atmosphere. In the case of a slowly varying refractive index, an approximation on the square root operator is necessary to end up with the phase screen method. We split the operator by means of the approximation proposed by Feit and Fleck [START_REF] Feit | Light propagation in graded-index optical fibers[END_REF] that provides a wide-angle approximation. It is given by

√ A + B + 1 ≈ √ 1 + A + √ 1 + B -1. (17)
This is the basic of the split-step algorithm where the freespace propagator √ 1 + A is computed in the spectral domain while the effect of refraction √ 1 + B -1 corresponds to a multiplication in the spatial domain by a phase screen. A higher-order approximation valid at wider angles can be applied [START_REF] Lin | A higher-order split-step Fourier parabolicequation sound propagation solution scheme[END_REF].

B. Numerical scheme

To simulate the propagation, the domain r ∈ [r 0 , r max ] is discretized in N r points, such that r = r 0 + p r ∆r for p r ∈ {1, ..., N r }. As with other split-step methods [START_REF] Levy | Parabolic Equation Methods for Electromagnetic Wave Propagation[END_REF], the computation is performed marching on in distance. Between two consecutive cylinders, we propagate the potential through a homogeneous medium using the spectral propagator. Then we apply a phase screen to account for the refraction index variations.

At each iteration p r , the propagation, the apodization and the local refractive index are considered using the 5 steps described below:

1) The potential Ψ on the cylinder at distance r 0 + p r ∆r is expressed in the spectral domain by means of the transform U defined by (8).

2) The spectrum is multiplied by the spectral operator Λ pr,q θ ,qz representing the propagation from p r ∆r to (p r + 1)∆r given by Λ pr,q θ ,qz = H

κ (k r (r + ∆r)) H

(2)

κ (k r r) r + ∆r r , (18)
with r = r 0 + p r ∆r.

3) The spectrum is expressed in the spatial domain with the transform U -1 . 4) To remove reflections over the top boundary, an apodization is applied with a Hanning window on the top half area in the spatial domain. This amounts to a term by term multiplication represented by a diagonal operator W . 5) A phase screen is applied which corresponds to a multiplication by R pr,p θ ,pz = exp (jk 0 (n pr,p θ ,pz -1)∆r) ,

where n pr,p θ ,pz is the refractive index at the position (p r ∆r, p θ ∆θ, p z ∆z). The refraction term can be seen as a perturbation of the homogeneous case.

To conclude, the propagation from p r ∆r to (p r + 1)∆r is simulated step by step as

Ψ pr+1 = RW U -1 ΛU Ψ pr (20)
where R, W , and Λ are diagonal operators, i.e. term-by-term multiplications. The final computation complexity of this method is of order

N r N θ N z log 2 N θ N z .
In the final step, the field is calculated from the potentials. To do so, the field spectrum is obtained from the potential spectrum by applying (2) in the spectral domain. Indeed, derivatives are easier to compute in the spectral domain. Finally, U -1 is applied to obtain the fields.

IV. SECTORAL PROPAGATION

When computing the propagation in a complex environment, two situations may occur. In the first one, you are only interested in the field within a limited angular width. In such a case, the most suitable solution in terms of computation time is to limit the computation domain to this angular sector and to impose absorbing boundary conditions [START_REF] Jensen | Modern Acoustics and Signal Processing[END_REF]. In the second one, you want to compute the complete field, without truncating any energy. This is our aim here. Rather than simulating the propagation in the entire space with θ ∈ [0, 2π[, the simulation can be limited to the sector where the power is located, denoted as the "propagation sector" thereafter. The azimuthal sector is chosen so that no energy reaches its boundaries. Doing so, any kind of boundaries can be imposed. Periodic boundary conditions are used because they maintain the validity of the spectral Fourier representation, which renders its computation achievable by FFT.

A propagation sector of width θ s in azimuth is chosen. The calculation domain is represented in Fig 2 . For the splitstep algorithm to remain unchanged, we impose θ s = 2π

Ns , where N s is an integer. Finally, the computation domain is

θ ∈ [-θ s /2, θ s /2[ and z ∈]0, z max [.
If periodic boundary conditions are imposed at ± θs 2 , the sectoral propagation is similar to the 2π-propagation except for a decrease of the angular period by a factor of N s . Thus, the truncation on θ amounts to a decimation of factor N s in the spectral domain. This means a N s times replication in the spatial domain as illustrated in Fig 3 . To avoid signal overlap, the sector should be chosen such that no power goes through its boundaries. A numerical criteria is proposed in the test sections.

Finally, only a slight change in the previous formulation is necessary. The expression of the propagator ( 18) is kept, except for the order of the Hankel function κ that becomes 2 ∆θ sin πNsq θ N θ for q θ = {0, ..., N θ -1}.

V. NUMERICAL TESTS OF THE 3D PROPAGATION

The aims of the simulations are to test the sectoral propagation and to show that the proposed 3D method achieves to model vertical and azimuthal effects. On the contrary, N×2D methods cannot take the lateral effects into account.

In the simulations, the choice of the steps ∆z and ∆θ depends on the Nyquist sampling theorem. The vertical and azimuthal samplings should be as [START_REF] Lin | Three-dimensional sound propagation models using the parabolic-equation approximation and the split-step Fourier method[END_REF] ∆z ≤ λ 2 sin(α max ) ,

and ∆θ ≤ λ 2r max sin(α max ) ,

where α max is the maximum divergence angle of the beam with respect to the propagation direction.

In the numerical simulations, the aim is to test the validity of the proposed 3D propagation method. For the sake of simplicity, we use a fixed azimuthal grid, even if this yields a large number of points in azimuth to satisfy the Nyquist condition at the final step, which are not necessary for the previous steps. Several schemes could be used to adapt the number of points in azimuth with range [START_REF] Lin | Three-dimensional sound propagation models using the parabolic-equation approximation and the split-step Fourier method[END_REF] [START_REF] Austin | the use of tessellation in three dimensional parabolic equation modeling[END_REF]. These schemes could be applied to the proposed formulation to reduce the computation time.

A. Validation of the sectoral propagation in free space

In order to test the sectoral propagation, a complex point source [START_REF] Deschamps | Gaussian beam as a bundle of complex rays[END_REF] is propagated and compared to its analytic expression. The frequency is 3 GHz. The complex source is located at Because a complex source corresponds to a Gaussian beam upon assuming the paraxial approximation, an angular sector of width 2 arctan 2 k0W0 contains 88.47 % of the beam power in the far field zone [START_REF] Goldsmith | Quasioptical Systems: Gaussian Beam Quasioptical Propagation and Applications[END_REF]. To guarantee that no power reaches the boundaries, the sector is chosen 5 times larger. This yields θ s = π 10 and N s = 20. The electric field at distance r max computed with the sectoral propagation is shown in Fig 4(a). The difference between the sectoral propagation and the analytic solution is shown in Fig 4(b). The maximum difference is of order -73dB. This difference is negligible. Therefore, the 3D sectoral propagation method is considered as successfully tested.

x s =
This simulation is performed on a standard desktop computer. The computation time of the sectoral propagation method is about 1.1 hours, whereas the total 2π propagation spends about 21 hours. The computation time is thus significantly reduced.

B. Test of a 45 • -slanted linear refractivity index model

The atmosphere effects are considered through the refractive index n and the modified refractivity M given by [START_REF] Pekeris | Accuracy of the earth-flattening approximation in the theory of microwave propagation[END_REF]

M = (n -1) × 10 6 + z R E , ( 23
)
where z is the height above the Earth's surface and R E is the Earth's radius.

The method is now applied on a 3D scenario with an inhomogeneous atmosphere characterized by a 45 • -slanted linear refractivity index.

The configuration parameters are as follows: the frequency is 3 GHz. The complex source is at an altitude of 500 m with W 0 = 3 m. The modified refractivity M varies linearly along the direction v = 1 √ 2 (ŷ + ẑ). The modified refractivity gradient along v is c = dM/dv = 1 M-unit/m. This value is chosen strong in order to obtain significant refraction effects, even at few kilometers.

The parameters of the simulation are as follows: r 0 = 2 km, r max = 12 km, ∆r = 500 m, z max = 1000 m, ∆z = 0.2 m, and N θ = 30000.

The initial and final fields are plotted in Fig. 5. The beam center of the final field is shifted in both altitude and azimuth due to refraction effects. The shifts of the beam center in z and θ are 35.81 m and 34.80 m respectively. According to ray theory [START_REF] Kravtsov | Geometrical Optics of Inhomogeneous Media[END_REF], this shift is given by

l shift = (r max -r 0 ) 2 c × 10 -6 2 . ( 24
)
The numerical application yields l shift = 35.36 m.

The ray approximation and the simulated values in the 2 directions are the same within an acceptable error. Therefore, the effects of the inhomogeneous atmosphere on vertical and horizontal directions are properly simulated. Note that the displacement along θ could not be simulated by a N × 2D method.

C. Tests on a 2D atmospheric duct

A more complex scenario with 3D effects along both vertical and azimuthal directions is presented. We consider a refractive index model identical in both directions with 2 ducts along the vertical and azimuthal directions. Therefore, similar effects are expected along both directions.

The refraction index in the z -θ plane is obtained by M = √ M z M θ , where M z and M θ are functions associated with the modified refractivity index along z and θ.

In the vertical direction, we consider a surface-based duct, which is modelled by a trilinear function M z , as illustrated in The plots match within 0.05 dB. Therefore, the differences of the fields along the 2 directions are negligible. The effects due to the ducts are the same in both directions. The method is successfully tested since 3D complex effects are simulated.

For a comparison, the result of a N×2D method is shown in Fig 11 . Only the vertical effects are accounted in this case.

VI. CONCLUSION

The propagation method introduced in this work is based on a discrete spectral representation of the propagation equation related to both vertical and azimuthal variables in cylindrical coordinates.

Firstly, the discretised form of the 3D formulation has been presented. A homogeneous atmosphere has been considered to express the problem by means of Hertz potentials. Then, the transform from the field to the spectrum has been explicitly given in discrete coordinates. A discrete spectral propagator has been derived. Secondly, an explicit numerical scheme for considering a slowly varying refractivity has been introduced. The phase screen method valid at wide-angles is applied in the spatial domain. The propagation is considered with Hankel functions. An apodization is finally applied on the top half area to remove reflections over the top boundary.

In many cases, the source is directional. Therefore, a sectoral propagation method has been introduced. The propagation is modified to fit with the reduced computation domain. Depending on the angular size of the propagation sector, the computation time can be drastically reduced.

Finally, three numerical tests have been presented. Firstly, the sectoral propagation and the 3D formulation have been validated by a comparison with an analytical solution. Then the method has been applied in a 3D scenario with an inhomogeneous atmosphere characterized by a 45 • -slanted linear refractivity index. The effects of the inhomogeneous atmosphere along vertical and horizontal directions are properly simulated. In the third simulation, a scenario with complex 3D effects along both vertical and azimuthal directions has been introduced. We have considered refractivity conditions with 2 atmospheric ducts along both the vertical and azimuthal directions. Note that the refractive effects have been successfully modelled along both directions. This method takes into account the azimuthal effects which is an advantage over N × 2D models.

Numerical tests have proven that the simulation time is significantly reduced when using the sectoral propagation. This is a first step but other speed-up approaches are currently investigated. Besides, the relief and depolarization effects should be introduced. Also, confrontations with measurements in the presence of azimuthally inhomogeneous atmosphere could be considered.

Fig. 1 .

 1 Fig. 1. Grids along z and θ on the initial cylinder at the distance r 0 and on a cylinder at distance r > r 0 .

Fig. 2 .

 2 Fig. 2. Propagation sector of azimuthal width θs = 2πNs (in pink).

Fig. 3 .

 3 Fig. 3. Decimation of a factor Ns in the spectral domain amounts to a Ns times duplications in the spatial domain (Ns = 8 in this figure).

2 ,

 2 y s = 0 m, z s = 1000 m, with W 0 = 1 m. The simulation parameters are: r 0 = 2 km, r max = 12 km, ∆r = 200 m, z max = 2000 m, ∆z = 0.2 m, and N θ = 1000. The propagation is in free space.

Fig. 4 .

 4 Fig. 4. Comparison between the 3D sectoral propagation and the analytic solution.

Fig 6 .

 6 The parameters are: M 0 = 330 M-units, z b = 950 m, z t = 100 m, z max = 2000 m, gradients c 0 = 0.118 M-units/m, c 2 = -1.0 M-units/m. To remove the discontinuities of the gradient that are not realistic, we smoothen the duct by using the moving averages method on 200 points. The resulting modified refractivity M z is shown in Fig 7(a). In the azimuthal direction, we use the same model mapped into cylindrical coordinates, as shown in Fig 7(b). Such large index variation in azimuth are not realistic in the troposphere, but are used to compare the 3D effect in θ and z. Finally, the modified refractivity index on the z -θ plane is obtained and the gradient of the refractivity index are shown in Fig 8 on one cylinder. For the simulation, the frequency is 3 GHz. The complex source is at an altitude of 1000 m with W 0 = 1 m. The simulation parameters are: r 0 = 2 km, r max = 12 km, ∆r = 500 m, z max = 2000 m, ∆z = 0.2 m, and N θ = 30000. Similar results along the 2 directions are expected since the ducts are the same along the 2 axes.

Fig. 5 .

 5 Fig. 5. Initial and final normalised electric fields (dB) about the direction of propagation.

Fig. 6 .Fig. 7 .Fig. 8 .

 678 Fig. 6. Surface based duct trilinear model.

Fig. 9 .Fig. 10 .

 910 Fig. 9. Final field of the simulation using the proposed 3D method.

Fig. 11 .

 11 Fig. 11. Final field of the simulation using N×2D method.