M G De

A M G De Lima

A R Da Silva

D A Rade

N Bouhaddi

Component mode synthesis combining robust enriched Ritz approach for viscoelastically damped structures

Keywords:

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

Introduction

The use of viscoelastic materials has been regarded as a convenient strategy in many types of industrial applications, where these materials can be applied either as discrete devices or surface treatments at a relatively low cost [START_REF] Nashif | Vibration damping[END_REF][START_REF] Rao | Recent applications of viscoelastic damping for noise control in automobiles and commercial airplanes[END_REF][START_REF] Samali | Use of viscoelastic dampers in reducing wind and earthquake induced motion of building structures[END_REF]. However, it presents some inherent drawbacks such as the influence of operational and environmental factors (frequency, temperature, pre-loads, etc.) [START_REF] Nashif | Vibration damping[END_REF]. Also, viscoelastic dampers (specially, surface treatments) are prone to induce considerable mass additions. A typical example is the application of viscoelastic constrained layer patches in compressor systems. In fact, having selected a viscoelastic material with optimum damping capability over the desired temperature and frequency range of compressor's operation, higher damping is achieved by increasing the viscoelastic treatment that results in an augmentation of the mass of the compressor [START_REF] De Lima | Vibration attenuation in refrigerators using constrained viscoelastic layers: preliminary finite element modeling and experimental assessment[END_REF]. This last feature leads to the necessity of performing optimization aimed at achieving the desired performance while complying with design and construction constraints. However, for finite element (FE) models of engineering structures incorporating viscoelastic materials composed by many thousands of degrees-of-freedom (DOFs) (for example an aerospace structure can be composed by more than 10 6 DOFs) the time to compute the exact evaluations during the iterative processes, performed on the full FE matrices, can become prohibitive [START_REF] De Lima | An efficient modeling methodology of structural systems containing viscoelastic dampers based on frequency response function substructuring[END_REF][START_REF] Balmès | Tools for viscoelastic damping treatment design: Application to an automotive floor panel[END_REF]. Those difficulties motivate the study reported herein with the intent to propose a general strategy to reduce the size of the viscoelastic structures for mainly several reasons: the cost of computation when one wants to extract the eigensolutions or to predict the behavior of the full structure, the optimization or model updating procedures and reliability-based optimization which require fast iterative techniques.

The difficulty of any model reduction procedure of viscoelastic systems lies on the fact that the viscoelastic stiffness matrix depends on frequency and temperature. Moreover, supplementary difficulty is to complete the representation basis in order to reduce truncation effects. One strategy is to generate a set of Ritz vectors capable of representing with precision the structural behavior under a wide variety of structural modifications. For example, Balmès and Germès [START_REF] Balmès | Tools for viscoelastic damping treatment design: Application to an automotive floor panel[END_REF] studied the possibility of using a constant Ritz basis to create parametric families of reduced models. Bouazzouni et al. [START_REF] Bouazzouni | Selecting a Ritz basis for reanalysis of frequency response functions of modified structures[END_REF] developed an optimal method to construct additional vectors by using the dynamic behavior of the structure before modification. The disadvantage of such approaches lies in the fact that a basis of reduction composed by a great number of residual static vectors is obtained, augmenting the computation effort involved in the condensation.

The substructuring or also known as Component Mode Synthesis (CMS) method currently plays a considerable role in the analysis of the viscoelastic systems, where a geometrically complex structure can be subdivided into components called substructures or superelements, which are analyzed and condensed separately, while preserving the junction's DOFs between the substructures. The substructures are represented by their modes (in the sense of Ritz vectors), including the normal, the rigid body, the static, and the interface modes. The structural modifications of a substructure make it possible to reanalyze the full system with a low cost of computation and make it possible to prepare each component independently of the others. This advantage enables the analysis of large-scale viscoelastic problems with the current data-processing resources.

Several CMS methods have been developed for forty years, and all of them make use of the vibration normal modes of the substructures. Depending on the boundary conditions applied to the substructure interfaces when these normal modes are obtained, the CMS methods can be classified into four groups: fixed interface methods [START_REF] Jr | Coupling of substructures for dynamic analysis[END_REF][START_REF] Hurty | Dynamic analysis of structural systems using component modes[END_REF][START_REF] Macneal | A hybrid method of component mode synthesis[END_REF][START_REF] Rubin | Improved component mode representation for structural dynamic analysis[END_REF]; free interface methods [START_REF] Macneal | A hybrid method of component mode synthesis[END_REF][START_REF] Rubin | Improved component mode representation for structural dynamic analysis[END_REF]; hybrid interface and loaded interface methods [START_REF] Rubin | Improved component mode representation for structural dynamic analysis[END_REF]. The variants of each group differ mainly in the choice of the supplementary Ritz vectors, and in the coupling procedure. In this paper one proposes an improvement of the classic Craig-Bampton Transformation (CBT) adapted for viscoelastic systems, by the introduction of residual static vectors associated to the external loads and the forces due to the viscoelastic damping.

In the remainder, after the various aspects related to the theoretical foundations, the description of a numerical application composed by a stiffened panel treated by a passive constraining damping layer demonstrates the effectiveness of the robust condensation strategy of viscoelastic systems.

Finite element models incorporating viscoelastic damping

Consider the following finite element equations of motion of a viscoelastic structure in the frequency domain containing N DOFs [START_REF] De Lima | An efficient modeling methodology of structural systems containing viscoelastic dampers based on frequency response function substructuring[END_REF]:

K e + K v (ω, T ) -ω 2 M Q (ω, T ) = F (ω) (1a) F (ω) = b u(ω); y(ω, T ) = cQ (ω, T ) (1b)
where M ∈ R N×N is the mass (symmetric, positive-definite) matrix, and K e ∈ R N×N and K v (ω, T ) ∈ R N×N are the stiffness matrices (symmetric, nonnegative-definite) corresponding to the purely elastic and viscoelastic substructures, respectively. Q (ω, T ) ∈ R N and F (ω) ∈ R N are, respectively, the vectors of displacements and external loads. y(ω, T ) ∈ R c is the vector of complex responses and u(ω) ∈ R f is the reduced vector of external loads. Matrices b ∈ R N×f and c ∈ R c×N are Boolean matrices which enable to select, among the d.o.f.s, those in which the excitation forces are applied and responses are computed, respectively. In expressions [START_REF] Nashif | Vibration damping[END_REF], the other forms of damping and the dissipation of the original structure are neglected.

It is widely known that the dynamic behavior of viscoelastic materials depend on a number of factors, among which the most relevant are the excitation frequency and the temperature [START_REF] Nashif | Vibration damping[END_REF]. Various mathematical models have been developed to represent this behavior and have been shown to be particularly suitable to be used in combination with FE discretization [START_REF] Bagley | A generalized derivative model for an elastomer damper[END_REF][START_REF] Golla | Dynamics of viscoelastic structures -A time domain finite element formulation[END_REF][START_REF] Mctavish | Modeling of linear viscoelastic space structures[END_REF][START_REF] Lesieutre | Time domain modeling of linear viscoelasticity using anelastic displacement fields[END_REF]. In this paper, as the interest is confined to frequency-domain analyses, the so-named Complex Modulus is used in combination with the Frequency-Temperature Correspondence Principle and the Elastic-Viscoelastic Correspondence Principle. According to the Complex Modulus approach, the dynamic behavior of viscoelastic materials in the frequency domain can be represented by introducing frequency-and temperature-dependent complex material moduli as follows [START_REF] Nashif | Vibration damping[END_REF]:

G(ω, T ) = G (ω, T ) [1 + iη G (ω, T )] (2)
where ω and T denote, respectively, the excitation frequency and the temperature of the viscoelastic material, G (ω, T ), η G (ω, T ) = G (ω, T )/G (ω, T ) and G (ω, T ) designate, respectively, the sonamed storage modulus, loss factor and loss modulus. It should be noted that definition (2) apply to both longitudinal and transverse moduli E(ω, T ) and G(ω, T ). However, it has been assumed by most authors [START_REF] De Lima | An efficient modeling methodology of structural systems containing viscoelastic dampers based on frequency response function substructuring[END_REF][START_REF] Balmès | Tools for viscoelastic damping treatment design: Application to an automotive floor panel[END_REF][START_REF] Bagley | A generalized derivative model for an elastomer damper[END_REF][START_REF] Golla | Dynamics of viscoelastic structures -A time domain finite element formulation[END_REF][START_REF] Mctavish | Modeling of linear viscoelastic space structures[END_REF][START_REF] Lesieutre | Time domain modeling of linear viscoelasticity using anelastic displacement fields[END_REF], that the Poisson ratio is independent from frequency and temperature in such a way that E(ω, T ) and G(ω, T ) are related to each other through the relation G(ω, T ) = E(ω, T )/ [2(1 + ν)]. Such an assumption, which has been adopted in this study, has been argued by [START_REF] Moreau | Identification de propriétés viscoélastiques de matériaux polymères par mesures de réponses en fréquences de structures[END_REF], who verified experimentally the variations of the Poisson ratio for PVC specimens.

The Frequency-Temperature Superposition Principle (FTSP), also known as Williams, Landell and Ferry (WLF) Principle [START_REF] Nashif | Vibration damping[END_REF], establishes an equivalence between the effects of the excitation frequency and temperature on the properties of a broad class of thermorheologically simple viscoelastic materials. This implies that the viscoelastic characteristics at different temperatures can be related to each other by changes (or shifts) in the actual values of the excitation frequency. This leads to the concepts of shift factor and reduced frequency, symbolically expressed as:

G(ω, T ) = G(ω r , T 0 ) = G(α T ω, T 0 ), η(ω, T ) = η(ω r , T 0 ) = η(α T ω, T 0 ) (3)
where T indicates an arbitrary value of the temperature, T 0 is a reference value of temperature, ω r = α T (T )ω is the reduced frequency, ω is the actual excitation frequency, and α T (T ) is the shift function.

Fig. 1 illustrates the FTSP, showing that having the modulus and loss factor of an arbitrary viscoelastic material for different temperature values, T -1 , T 0 , T 1 , if horizontal shifts along the frequency axis are applied to each of these curves, all of them can be combined into a single one, named master curve. The horizontal shift is given by α T and depends on the temperature.

Functions G(ω r ) and α T (T ) can be obtained from experimental tests for specific viscoelastic materials [START_REF] Nashif | Vibration damping[END_REF]. Drake and Soovere [START_REF] Drake | A design guide for damping of aerospace structures[END_REF] suggest analytical expressions for the complex modulus and shift factor for various commercial viscoelastic materials. Eqs. (4) represent the complex modulus and shift factor as functions of temperature and reduced frequency in the intervals 210 ≤ T ≤ 360 K and 1.0 ≤ ω ≤ 1.0 × 10 6 Hz, for the 3M TM ISD112 viscoelastic material, as provided by those authors. The 3M TM ISD112 is a rubber-like viscoelastic polymer which is provided by the manufacturer in the form of adhesive tapes [18].

G(ω r ) = 430 700 + 1200 × 10 6 1 + 3.241×(iω r /1543 000) -0.18 +(iω r /1543 000) -0.6847 (4a)

log(α T ) = -3758.4 × 1 T - 1 290 + 225.06 × log T 290 + 0.23273 × (T -290) (4b)
where T is the temperature in Kelvin, ω r is the reduced frequency in rad/s, and G is the complex modulus in N/m 2 . Fig. 2 depicts the standardized curves representing the variations of the storage and loss moduli, the loss factor and the shift factor as functions of the reduced frequency and temperature, respectively, as obtained from Eqs. (4a) and (4b).

According to the Elastic-Viscoelastic Correspondence Principle [START_REF] Christensen | Theory of viscoelasticity: An introduction[END_REF] the derivation of the FE model accounting for the viscoelastic behavior can be carried out in two distinct phases: first, the element and global stiffness matrices are obtained by considering pure elastic behavior (hence, frequency-and temperatureindependent material moduli), accounting for the strain state assumed; then, the material moduli are modified to account for the viscoelastic behavior (according to the model expressed by Eq. ( 2), in the present case). This implies that G(ω, T ) can be factored-out of the matrix K v (ω, T ) = G(ω, T ) Kv (where Kv is a constant ma- trix), which can be combined with Eqs. [START_REF] Nashif | Vibration damping[END_REF] to produce the receptance or frequency response function (FRF) matrix: where Z (ω, T ) = K e + G(ω, T ) Kvω 2 M is the so-named dynamic stiffness matrix.

H (ω, T ) = cZ (ω, T ) -1 b (5)
The difficulty in predicting the vibration response in the time domain and performing eigenvalue analysis for the viscoelastic systems comes from the fact that the viscoelastic stiffness matrix depends on frequency and temperature. As a result, one has a nonlinear eigenvalue problem that must be solved iteratively [START_REF] Daya | A numerical method for nonlinear eigenvalue problems: Application to vibrations of viscoelastic structures[END_REF]. Some procedures for dealing with this problem have been suggested, such as the recent contributions by Palmeri and Ricciardelli [START_REF] Palmeri | Fatigue analyses of buildings with Viscoelastic dampers[END_REF], in which the eigenvalue problem of viscoelastic systems has been derived in the time domain by the concept of modal relaxation functions. Yuan and Agrawal [START_REF] Yuan | A numerical scheme for dynamic systems containing fractional derivatives[END_REF] and by Wagner and Adhikari [START_REF] Wagner | Symmetric state-space method for a class of nonviscously damped systems[END_REF] proposed an alternative state-space approach for the time-domain analysis of viscoelastic systems. Other alternatives have been proposed based on the adoption of particular representations for the frequency-dependent behavior of the viscoelastic materials, such as the Fractional Derivative [START_REF] Bagley | A generalized derivative model for an elastomer damper[END_REF], Golla-Hughes-McTavish [START_REF] Golla | Dynamics of viscoelastic structures -A time domain finite element formulation[END_REF][START_REF] Mctavish | Modeling of linear viscoelastic space structures[END_REF] and Anelastic Displacement Field [START_REF] Lesieutre | Time domain modeling of linear viscoelasticity using anelastic displacement fields[END_REF] models, which enable to transform the equations of motion of the viscoelastic systems in the time-domain into statespace forms, with frequency-independent state matrices, at the expense of a typically high increase in the order of the system matrices.

The interest herein is focused on frequency-domain responses, and with this aim, Eq. ( 5) can, in principle, be directly used for calculating the steady-state harmonic responses of the viscoelastic structures in the frequency domain. However, such a procedure can be unfeasible when large-scale finite element models are dealt with. In this case, the computations can be alleviated by using model reduction techniques.

Robust condensation of viscoelastic systems for design procedure

In the case of complex engineering structures of industrial interest, FE models are usually constituted by a large number of DOFs (hundreds of thousand or even millions). In such cases, it becomes practically impossible to compute the FRFs directly from Eq. ( 5), owing to the prohibitive computation times and storage memory required. This fact motivates the use of model reduction procedures, which aims at reducing the model dimensions (and the associated computational burden), while keeping a reasonable predictive capacity of the numerical models. This can be done based on the assumption that the exact responses, given by the resolution of Eqs. [START_REF] Nashif | Vibration damping[END_REF], can be approached by projections on a reduced vector basis as follows:

Q (ω, T ) = T Q (ω, T ) (6)
where T ∈ C N×NR is the transformation matrix formed columnwise by a vector basis, Q (ω, T ) ∈ C NR are generalized coordinates, and NR N is the number of reduced vectors in the basis. The generalized coordinates representing the contribution of each column of T are chosen arbitrarily in which the reduce model provides a reasonable predictive capacity into a frequency bandwidth. Frequency band of interest is taken into account by computing a number of normal modes and retaining those below a certain frequency (1.5 times the last frequency of interest is typically).

By considering Eqs. ( 1) and ( 6), the transfer function ( 5) can be rewritten as:

Ĥ(ω, T ) = c Ẑ(ω, T ) -1 b (7) where Ẑ(ω, T ) = T T K e T + G(ω, T )T T Kv T -ω 2 T T MT . Ẑ(ω, T ) ∈
R NR×NR is the reduced dynamic stiffness matrix that must be computed frequency by frequency and inverted in a direct way by using efficient numerical algorithms. For models containing viscous or structural damping, it is relatively common to use constant projection basis formed by the eigenvectors of the associated conservative structure, as the mass and stiffness matrices are invariant [START_REF] Nashif | Vibration damping[END_REF]. However, for viscoelastic systems, the selection of the reduction basis is more delicate as this condition does not hold. Owing to the dependence of the stiffness matrix with respect to frequency and temperature, the reduction basis should be able to represent the changes of the dynamic behavior as frequency and temperature vary. Three procedures have been adopted regarding the computation of the reduction basis: (i) one can simply neglect this dependence by considering the stiffness matrix as being constant [START_REF] Balmès | Tools for viscoelastic damping treatment design: Application to an automotive floor panel[END_REF]. In this case, the reduction basis is also constant; (ii) one can use a constant reduction basis obtained by the resolution of the nonlinear eigenvalue problem associated to a frequency-and temperature-dependent stiffness matrix [START_REF] Daya | A numerical method for nonlinear eigenvalue problems: Application to vibrations of viscoelastic structures[END_REF]; (iii) one can use an iterative method for the re-actualization of the reduction basis according to frequency [START_REF] Balmès | Tools for viscoelastic damping treatment design: Application to an automotive floor panel[END_REF].

In this work, the strategy proposed consists in using a reduction basis formed by a constant modal basis (named herein nominal reduction basis), enriched by static residual vectors or equivalently static responses to account for the static effects of the modal truncation. These static responses are computed by using the tangent stiffness matrix representing the static behavior of the viscoelastic materials. As can be seen in Fig. 2(a), in the low frequency range, as the modulus and loss factor curves are prolonged by asymptotes, the extrapolation leads to the real values G 0 and η 0 = 0. On the other hand, for high frequencies, the extrapolation gives the complex values G ∞ and η ∞ . The tangent stiffness matrix can thus be calculated as follows [START_REF] Balmès | Tools for viscoelastic damping treatment design: Application to an automotive floor panel[END_REF]:

K 0 = K e + G 0 Kv . (8)
The nominal basis can be obtained by the resolution of the eigenvalue problem:

(K 0 -λ i M ) φ i = 0 i = 1, . . . , N φ 0 = φ 1 φ 2 . . . φ NR , Λ 0 = diag (λ 1 , . . . , λ NR ) . (9)
This basis is enriched by introducing the residues formed by the static displacements associated to external forces,

R = K -1 0 b (10)
that must be further completed by the residual vectors associated to the viscoelastic damping forces:

R 0 v = K -1 0 Kv φ 0 . (11)
As detailed in Reference [START_REF] Balmès | Tools for viscoelastic damping treatment design: Application to an automotive floor panel[END_REF], these residuals are interpreted as the columns of the flexibility matrix of the associated undamped system, associated to the coordinates of application of two types of forces to it: the external excitation forces and the damping forces. These latter can be better understood by examining [START_REF] Nashif | Vibration damping[END_REF], noting that the term involving the viscoelastic behavior can be moved to the right-hand side, where it plays the role of additional forces applied to the associated conservative structure. Thus, the enriched basis of reduction for the viscoelastic system is given as follows:

T 0 = φ 0 R R 0 v . (12)
It is important to note that the time required to compute the exact matrix K -1 0 for FE models composed by a large number of DOFs can become prohibitive. The numerical technique considered to be more accurate and time-effective consists in perform the Cholesky decomposition into the product of a lower triangular matrix and its conjugate transpose easier to be inverted [START_REF] Lay | Linear algebra and its applications[END_REF].

Experience has demonstrated that the nominal basis [START_REF] Bagley | A generalized derivative model for an elastomer damper[END_REF] can be used to reduce the viscoelastic systems with reasonable accuracy, but is not capable of representing the modifications of the dynamic behavior provoked by the modifications which must be introduced into the model during iterative optimization or model updating processes. This means that the basis [START_REF] Bagley | A generalized derivative model for an elastomer damper[END_REF] should, in principle, be updated successively to guarantee a satisfactory accuracy of viscoelastic model reduction while the model evolves during iterations. However, this process would involve costly computations. To cope with this difficulty, the strategy suggested herein consists in performing a further enrichment of the nominal basis T 0 by a set of residual vectors calculated based on the knowledge of the parametric modifications of the viscoelastic zones.

Robust static residual vectors

The robust condensation strategy lies in its use of a priori information of the nature and localization of the potential design modifications, while keeping the amplitude of the modification as a variable [START_REF] Masson | Component mode synthesis (CMS) based on an enriched Ritz approach for efficient structural optimization[END_REF]. In such a way, for viscoelastic systems, the objective is to construct a set of robust frequency-independent Ritz vectors, to complete the standard reduction basis [START_REF] Bagley | A generalized derivative model for an elastomer damper[END_REF]. By considering the Eqs. ( 1) and ( 5), for the modified structural configuration on the viscoelastic zones, the dynamic equilibrium equation in the frequency domain can be written as follows:

[ Z v (ω, T ) + Z (ω, T )] Q (ω, T ) = F (ω) (13)
where Z v (ω, T ) = K v (ω, T )-ω 2 M v is the variation of the dynamic stiffness matrix associated to the viscoelastic modifications. Eq. ( 13) can be interpreted as the dynamic equilibrium equation of the nominal viscoelastic model, subjected to the forces F (ω, T ) associated to the structural modifications:

Z (ω, T )Q (ω, T ) = F (ω) + F (ω, T ) (14a)
where:

F (ω, T ) = -Z v (ω, T )Q (ω, T ). (14b)
At this point, it is important to consider that, in the context of the present study, the modifications are introduced on a set of geometrical design parameters associated to the viscoelastic zones that intervene in a rather complicated nonlinear fashion in those dynamic matrices. In an attempt to compute the modified matrices involved in the further condensation, it becomes interesting to perform a parameterization of the FE model, which is understood as a means of making them to appear explicitly in the finite element mass and stiffness matrices. This procedure enables to account for such modifications and/or uncertainties in the values of the design parameters in a straightforward way during iterative processes. Also, it facilitates, to a large extent, the evaluation of the sensitivities of the responses with respect to the design parameters [START_REF] De Lima | Sensitivity analysis of frequency response functions of composite sandwich plates containing viscoelastic layers[END_REF][START_REF] De Lima | Sensitivity analysis of viscoelastic structures[END_REF][START_REF] Guedri | Robust design of viscoelastic structures based on stochastic finite element models[END_REF]. In general, the parameters of mass and stiffness intervene nonlinearly into the FE matrices, and after manipulations the perturbed matrices having the design parameters factored-out can be expressed symbolically as follows:

K v (ω, T ) = G(ω, T ) zones i=1 Kvi , M v = zones i=1 M vi (15a) Kvi = p i p i α α(p i ) α α Kvi , M vi = p i p i β β(p i ) β β M vi . (15b)
In the equations above, symbol indicates matrix assembling, subscript i indicate the modified zones (most frequently encompassing various viscoelastic elements) in which the parameter p i intervenes. Kvi and M vi designate the stiffness and mass matrices corresponding to those viscoelastic zones, respectively. These later are decomposed according to Eq. (15b) by identifying the matrices α Kvi and β M vi from which the αand β-order exponentials in the parameter p i can be factored-out.

As can be seen in Eq. (14b) the vector of forces associated to the structural modifications depends on the dynamic response of the modified structure, Q (ω, T ). Since this response is unknown a priori, those forces cannot be computed exactly. The essential concept of the robust condensation is expressed in the following [START_REF] Masson | Component mode synthesis (CMS) based on an enriched Ritz approach for efficient structural optimization[END_REF]: (i) Eq. (14b) is used to generate a vector of forces which, even though it does not contain the exact forces associated to the modifications, will at least represent a subspace containing these vectors. This is accomplished by introducing the response of the nominal system into Eq. (14b); (ii) the resulting vector of forces is used to generate static responses, once again on the basis of the nominal model; (iii) the first two steps are repeated for each design parameter subjected to modifications. Many types of responses may be introduced into Eq. (14b), including the normal modes and the sensitivity vectors. If the vector introduced is composed by a truncated basis of normal modes, Eq. (14b) assumes the following form:

F (ω, T ) ≈ -Z v (ω, T )φ 0 Q (ω, T ). ( 16
)
For example, for a parameter p i intervening in the viscoelastic mass and stiffness matrices, the basis of forces can be expressed under the following form:

F p i = F M vi p i F Kvi p i ( 17
)
where F M vi p i = M vi φ 0 Λ 0 and F Kvi p i = G 0 Kvi φ 0 are the basis of forces associated to the viscoelastic mass and stiffness modifications.

After obtaining the basis of forces, one can calculate a series of static responses of the nominal viscoelastic system based on the tangent stiffness matrix as follows:

R p i = K -1 0 F p i , i = 1, 2, . . . , np. ( 18
)
The final robust condensation basis taking into account a priori knowledge of the viscoelastic modifications can be expressed as follows:

T = T 0 R (19)
where R = R p 1 R p 2 . . . R p np , and np indicates the number of design parameters to be modified. The residue matrix R is not necessarily of maximum rank. Thus, with the aim of obtaining a limited number of independent residue vectors, it is appropriate to select the dominating directions of this basis, which can be done by performing the Singular Value Decomposition (SVD) of R to identify its dominant singular values.

Fig. 3 illustrates a comparison between a cycle of optimization and/or model updating processes by using a standard reduction and the proposed robust condensation strategy. The robust strategy is used to approximate the behavior of the modified viscoelastic structures without the re-actualization of the nominal basis of reduction, leading to a drastic reduction of the time required for computing the responses. Moreover, this approach can be advantageously adapted to several other structural domains based on iterative processes: stochastic structural dynamics [START_REF] Guedri | Robust design of viscoelastic structures based on stochastic finite element models[END_REF], nonlinear mechanics, and reliability-optimization-based design.

In the case of global modifications, for which the whole structure fully treated by viscoelastic materials is modified by using the same level of perturbation, the standard condensation procedure by using the reduction basis ( 12) is robust and do not need to be improved. In this case, there is only a frequency shift between the nominal model and the perturbed model, so that transformation [START_REF] Bagley | A generalized derivative model for an elastomer damper[END_REF] can represent correctly the perturbed viscoelastic model. But when the local modifications are introduced that is the case of structures partially treated by passive constraining layer damping, experience shows that the precision of this solution degenerates rapidly with the amplitude of the perturbation, demonstrating the interest in using the robust basis [START_REF] Christensen | Theory of viscoelasticity: An introduction[END_REF].

CMS combining robust condensation for viscoelastic systems

The objective of the substructuring methods is to calculate the dynamic behavior of a complex engineering structure, starting from the knowledge of the dynamic behavior of its substructures [START_REF] Jr | Coupling of substructures for dynamic analysis[END_REF][START_REF] Hurty | Dynamic analysis of structural systems using component modes[END_REF]. The principal differences between the substructuring methods are the choice of the basis of reduction for each substructuring, and the assembly technique between them. The principal sources of errors when implementing a substructuring method are: (i) the use of incomplete Ritz basis; (ii) errors in the discretization process of the geometry of an interface between two components; (iii) and the inequality of displacements and the balance of forces over an interface.

Fig. 4 illustrates two substructures finite element mesh whose DOFs can be partitioned in two subsets: an internal domain composed by all internal nodes belonging to the mesh (exception of the right edge), represented by the subscript i, and the nodes of the right edge defining the junction, represented by the subscript j.

The DOFs of the displacement vector Q (ω, T ) of a viscoelastic substructure SS2 are partitioned into two subsets: one related to the junction DOFS, Q j (ω, T ), and the others associated to the interior DOFs, Q i (ω, T ), according to the expression: The dynamic equilibrium equation of the substructure SS2 characterized by the tangent stiffness matrix (according to Eq. ( 8) in Section 3), can be written by taking into account the partition (20) as follows:

Q (ω, T ) = Q j (ω, T ) Q i (ω, T ) . (20)
K 0jj K 0ji K 0ij K 0ii -Λ M jj M ji M ij M ii Q j Q i = F j 0 ( 21
)
where Λ = diag(λ 1 , . . . , λ NS ), NS is the number of DOFs of the substructure.

The assembly properties of the reduced mass and stiffness matrices depend on the type of selected modes and the assembly process. In the case of assembly level of reduced matrices of several viscoelastic substructures, in this paper, one uses the method suggested by Craig and Bampton [START_REF] Jr | Coupling of substructures for dynamic analysis[END_REF], which imposes the continuity of displacements on the interface between the substructures, according to the following expression (for the system illustrated in Fig. 4):

Q (SS1) j = Q (SS2) j (ω, T ) . (22)
In dynamic analysis of complex viscoelastic structures with a large number of DOFs, the Craig-Bampton can be used to drastically reduce the overall number of DOFs when combined with component mode synthesis technique. The condensation basis known as traditional Craig-Bampton Transformation (CBT) is composed by the vectors associated to the forced modes enriched by normal modes for a substructure with blocked junction. The constraint modes correspond to the static response of the substructure when a unit displacement is imposed on one selected DOF of the contour, and the others DOFs are blocked. For a viscoelastic substructure, they form the solution of the following equation:

K 0jj K 0ji K 0ij K 0ii I jj G ij = F j 0 ( 23
)
where I jj represents a unit vector at DOF j, that is a vector with zero entries except at the DOF j that is equal to one. It represents a unit displacement imposed on the DOF j of the junction. Vector G ij represents the static displacements due to zero forces at the interior DOFs. The resolution of Eq. ( 23) yields:

G ij = -K -1 0ii K 0ij I jj , F j = K 0jj -K 0ji K -1 0ii K 0ij I jj . ( 24
)
The normal modes correspond to the eigensolutions of the viscoelastic substructure. There are three different types of normal modes according to the boundary conditions imposed on the interface [START_REF] Hurty | Dynamic analysis of structural systems using component modes[END_REF]. In this work, one uses the modes associated to the blocked interface, ψ i , obtained by the resolution of the homogeneous viscoelastic problem associated to the internal DOF i, for which the DOF j is subjected to the boundary conditions. The vector representing the generalized coordinates are then composed by the physical coordinates Q j and the modal coordinates Q m associated to the normal modes with blocked interfaces ψ. The CBT transformation for a viscoelastic substructure is defined as follows:

Q SS2 = Q j Q i = I jj 0 G ij ψ Q j Q m (25a)
where:

T CBT = I jj 0 G ij ψ . (25b)
The reduced problem is obtained by substituting Eq. (25b) into Eq. (1a) and premultiplying the resulted system equation by the transpose of the basis T CBT . However, for the same reason as explained in the previously section, for a viscoelastically damped substructure, the basis (25b) must be enriched by static residues associated to the external loads, and the forces associated to the viscoelastic modifications. The normal modes of the viscoelastic substructure chosen to construct the first order static residues are the modes associated to the blocked interface, and the residual vectors representing the static displacements related to the external forces, R, according to Eq. ( 10), and the residues associated to the viscoelastic modifications with blocked junction, R 0 v and R , calculated according to Eqs. ( 11) and (18), respectively. Thus, the robust enriched CBT, T CBTE , for a viscoelastic substructure is:

T CBTE = T 0 CBTE I jj 0 0 0 G ij ψ R R 0 v 0 R . (26)
The strategy for improving the robustness of the classical CBT basis (25b) with respect to structural viscoelastic modifications to obtain the basis [START_REF] De Lima | Sensitivity analysis of frequency response functions of composite sandwich plates containing viscoelastic layers[END_REF], is similar to that described in Section 3.1, in which the optimal set of Ritz vectors R are constructed based on the introduced modifications on the nominal viscoelastic model, to complete the standard enriched CBT basis T 0 CBTE for the viscoelastically damped substructure.

Numerical application

In this section, a numerical application is presented to illustrate the main features of the proposed robust condensation strategy of viscoelastic systems. also indicated on Fig. 5(a), is composed by 10 viscoelastic patches, each one comprising 16 three-layer sandwich plate elements developed according to References [START_REF] De Lima | Sensitivity analysis of viscoelastic structures[END_REF][START_REF] Guedri | Robust design of viscoelastic structures based on stochastic finite element models[END_REF]. The FE model of the partially viscoelastically treated panel contains 6840 DOFs. The geometric dimensions are: internal radius: 938 mm; length: 720 mm; arc length: 680 mm; thicknesses of the panel and the stringers: 1.5 mm and 0.75 mm, respectively; height of the stringers: 30 mm. The material properties for both panel and stringers are: Young modulus E = 2.1 × 10 11 N/m 2 ; mass density ρ = 7800 kg/m 3 ; Poisson's ratio, υ = 0.3. The material properties of the constraining layer are the same as those of the stiffened panel, and for the viscoelastic core, one uses the modulus function of the 3M TM ISD112 (ρ = 950 kg/m 3 ), as shown in Fig. 2(a). The nominal values of the thicknesses of the viscoelastic and constraining layers are 0.0254 mm and 0.5 mm, respectively.

Evaluation of the robust condensation basis

The interest here is to verify the robustness of the reduced viscoelastic model generated by the use of the proposed condensation method described in Section 3.1. The first test is intended to evaluate the nominal enriched basis of reduction by using the static residues associated to the external loads and the viscoelastic forces. The computations consisted in obtaining the driving point FRFs associated to point P indicated on Fig. 5 10)); T 03 = φ 0 R R 0 v (60 eigenvectors, one residual vector computed by Eq. ( 10), 54 residual vectors computed according to Eq. ( 11) after SVD filtering). The residues R 0 v were computed based on the largest singular values, for which the relation σ 1 /σ i ≤ 1 × 10 5 , for i = 1 to 60 holds. The interest now is to evaluate the robustness of the nominal basis further enriched to account for structural modifications introduced into the nominal model, according to Eq. ( 19). The modification considered consists in increasing the thickness of the constraining layer of the nominal system by 90%. Fig. 8 enables us to compare the amplitudes of the frequency response functions for the nominal and perturbed damped systems both computed using the same reference reduction basis, without further enrichment for structural modifications. It can be seen that the dynamic behavior of the perturbed viscoelastic damped system does not differ too much from that of the nominal system, although viscoelastic damping levels are strongly influenced.

In Fig. 9(a), the FRF of the perturbed reference system is compared to the counterpart computed by using the nominal basis T 03 , containing 115 vectors. The observed differences lead to conclude that this basis is not capable of accurately representing the changes of the dynamic behavior induced by the structural modifications introduced. agreement between the amplitudes of the dynamic responses of both systems. This leads to conclude that the reduction basis T 4 is robust enough to represent the dynamic response of the perturbed viscoelastic structure in the frequency band of interest, yielding a significantly more accurate prediction. Fig. 9(c) confirms that the use of residual vectors associated to the structural modifications is sufficient to represent with accuracy the dynamic behavior of the modified viscoelastic system.

Component mode synthesis

To verify the dynamic substructuring condensation strategy of viscoelastic structures, one uses the substructure definition represented in Fig. 5 By using the enriched CBT basis ( 26), the nominal model of the first substructure is condensed to a substructure having 221 vectors described as follows:

T 0 CBTE = I 180×180 0 180×20 0 180×1 0 180×20 G 4170×180 ψ 4170×20 R 4170×1 R 0
where G 4170×180 represents the constrained modes; ψ 4170×20 represents the normal modes associated to the blocked interface; R 4170×1 is the residues associated to the external loads; and R 0 v 4170×20 represents the residues associated to the viscoelastic forces.

In the same manner, the condensed model of the second substructure (purely elastic) is composed by 190 vectors, for which the classical CBT is expressed as:

T CBT = I 180×180 0 180×10 G 2340×180 ψ 2340×10 (190 vectors). (28)
Fig. 10(a) shows the FRF amplitudes of the reference system and those obtained for the condensed model by using the enriched CBT. It is shown that the enrichment of the classical CBT basis of reduction associated to the external loads and the viscoelastic effects improves the results in the frequency band of analysis.

The modifications imposed on the nominal model consist in perturbing the viscoelastic zones of the damped substructure, by increasing the constraining layer thickness by 90%. The reduction of the modified system is made by using the nominal basis of reduction associated the nominal system (compound by 221 vectors), enriched by static residues evaluated on the modified matrices, to take into account the modifications on the viscoelastic parameters.

T CBTE = T 0 CBTE R (246 vectors).
(29) Fig. 10(b) enables us to compare the dynamic responses of the perturbed systems computed for the reference and reduced systems, respectively. Again, by this figure, one can conclude that the first order static correction associated to the viscoelastic modifications is robust enough to represent correctly the dynamic response of the perturbed system in the frequency band of interest.

To provide a sense of the additional computation effort associated with the enrichment process, Table 1 provides the comparison between the direct and CMS reduction methods. The most immediate use of the enrichment process is the increase of the total time required to compute the reduced system response [START_REF] Bouazzouni | Selecting a Ritz basis for reanalysis of frequency response functions of modified structures[END_REF]. Also, the assembly of the matrices and their reduction take more time of the total computing time for the direct reduction, which motivates the use of a parameterization of the FE model, leading to significant time savings in iterative processes. The results shown in Table 1 demonstrate that the CMS method is a little bit more costly than a direct method in term of computation but definitely more economic than an exact reanalysis when dealing with very large models having a large number of subassemblies. Moreover, one can take advantage of this enriched substructuring scheme when only one component incorporating surface viscoelastic treatment is subject to local modifications. In this case, there is no need for re-assembling or projecting on the reduction basis [START_REF] De Lima | Sensitivity analysis of viscoelastic structures[END_REF] the matrices of the others elastic components. Thus, the iterative processes would give better results in terms of computing time.

Concluding remarks

A robust condensation procedure combined with a substructuring technique, intended to be used for dealing with large-scale viscoelastically damped systems was suggested and evaluated. The original aspects of the procedure lies in the adaptation of the concept of robust condensation, initially developed for undamped structures, for systems containing viscoelastic materials, and the extension of the proposed methodology to the case of CMS method. This approach allows design parameters of a superelement to be modified, for example in the context of an optimization and/or model updating processes, without the necessity to perform a complete superelement analysis at each point in parameter space. Thus, the improved model reduction transformation can be prepared in advance for each substructure and then used directly during the iteration processes to avoid the exact recalculation of the modified zones.

For viscoelastic structures the proposed robust condensation methodology can be integrated in a variety of component mode synthesis techniques and we have illustrated its use in the context of the Craig-Bampton condensation procedure. Moreover, this approach is particularly effective when dealing with very large-scale FE viscoelastic models having a large number of substructures. At the present time, proposed design process is not constrained solely to the viscoelastic elements, and can already be used to approximate reanalysis of linear elastodynamic behavior to substructure technology.

An academic example derived from the aeronautic engineering interest was used to illustrate the efficiency of the robust condensation procedure in representing the modified behavior of the perturbed viscoelastic zones. The obtained results demonstrated the effectiveness of the two robust condensation strategies mainly in terms of the drastic reduction of the number of DOFs, which demonstrates that the suggested technique is well adapted to be applied to complex viscoelastic structures.

Fig. 1 .

 1 Fig. 1. Illustration of the Frequency-Temperature Superposition Principle (FTSP).

Fig. 2 .

 2 Fig. 2. Master curves (-G; -• -η) for the 3M TM ISD112 viscoelastic material.

Fig. 3 .

 3 Fig. 3. Block-diagram of a reanalysis process combining robust condensation.

Fig. 4 .

 4 Fig. 4. Illustration of two substructures FE mesh with partitioned DOFs.

Fig. 5 Fig. 5 .

 55 Fig. 5. FE model of the test structure (a); Substructure definition (b).

Fig. 6 .

 6 Fig. 6. FRF amplitudes of the reference and reduced systems by using the nominal basis T 01 (a), T 02 (b) and T 03 (c).

  (a), assuming a temperature of the viscoelastic material of 25 • C, and a frequency band of interest [135-210 Hz]. To verify the direct condensation, one considers the following nominal basis: T 01 = φ 0 (60 eigenvectors); T 02 = φ 0 R (60 eigenvectors, plus one residual vector computed by Eq. (

Fig. 6 (

 6 Fig.6(a), (b) and (c) show the FRFs computed by using the three nominal basis, as compared to the amplitudes of the response computed by using a reference basis formed by a far larger

  Fig. 9(b) enables us to compare the FRF of the perturbed reference system to the counterpart computed by using the basis T 4 = T 03 R containing 146 vectors (including

Fig. 7 .

 7 Fig. 7. Comparison between the relative errors on the FRF amplitudes computed for the three reduction basis (a); relative errors computed for the basis T 03 (b).

Fig. 8 .

 8 Fig. 8. FRF amplitudes computed for the nominal and perturbed systems using the reference reduction basis.

  (b). The first substructure is composed by 551 shell finite elements, having 4350 DOFs (180 DOFs of junction), whose 160 finite elements are treated by 3M TM ISD112 viscoelastic material with the temperature of 25 • C, and the FE model of the second substructure is composed by a total of 377 finite elements, having 2520 DOFs, without viscoelastic treatment.

Fig. 9 .

 9 Fig. 9. FRF amplitudes computed for the nominal and perturbed systems using different reduction basis: (a) T 03 ; (b) T 4 .

Fig. 10 .

 10 Fig. 10. FRF amplitudes computed for the reference and reduced systems by using the enriched Craig-Bampton transformation.

Table 1

 1 Computing time for direct and CMS reduction methods.

		Direct method				CMS method
		T 01 (s)	T 02 (s)	T 03 (s)	T 4 (s)	T 0 CBTE (s)	T CBTE (s)
	Matrix assembly and reduction	18.3	18.4	19.9	3.3 a	14	3.8 a
	Time to compute the reduced system (7)	2.8 s	2.9	3.7	5.8	17	19.9
	Total computing time	21.1	21.3	23.6	9.1	31	23.7

a Time to compute only the residues associated to the structural modifications and to perform the reduction.

modifications). This time, one can observe a very satisfactory

Acknowledgements

The authors wish to thank Brazilian Research Council -CNPq for the continued support to their research work, especially through research projects 480785/2008-2 (A.M.G. de Lima) and 310524/2006-7 (D.A. Rade); Minas Gerais State Agency FAPEMIG, CAPES Foundation, and French agency COFECUB for the support to the international cooperation project with the Applied Mechanics Laboratory R. Chaléat, FEMTO-ST.


Sameer Alam
email: s.alam@adfa.edu.au

Supatcha Chaimatanan
email: supatcha@gistda.or.th

Daniel Delahaye
email: delahaye@recherche.enac.fr

Eric Feron
email: feron@gatech.edu

A Distributed Air Traffic Flow Management Model for European Functional Airspace Blocks

Functional airspace blocks (FAB) concept is adopted by the European airspace to allows cooperation between airspace users to manage the air traffic flow, while ensuring efficiency, safety, and fairness without constraints of geographical boundaries. This integration of airspace introduces more flexibility to manage aircraft trajectory and airspace usage. In this paper, we propose a distributed air traffic flow management model to address four-dimensional trajectory planning over the European FAB. The goal is to enable effective information sharing between airspace blocks in strategic planning, and to minimize interaction between trajectories. The proposed model and overall methodology is implemented and tested with a full day air traffic data over the European airspace. Interaction-free 4D trajectories are produced within computational time acceptable for the operational context, which shows the viability of the distributed model and interaction minimization approach for effective FAB implementation.

I. INTRODUCTION

The main objective of an Air Traffic Flow Management (ATFM) system is to ensure safety, reduce delays, and balance demand and capacity among different components of the air transportation system [10]. In European airspace, ATFM activities are carried out by Eurocontrol's Network Manager Operations Centre (NMOC, previously called CFMU) which acts as a coordinator between air navigation service providers (ANSPs) and key stake holders, such as airlines, airports, and military.

However, fragmentation of airspace along national boundaries creates structural inefficiencies, sub-optimal flight routing and makes implementation of ATFM strategies quite challenging [START_REF] Button | Single european sky and the functional airspace blocks: Will they improve economic efficiency[END_REF]. To address these challenges, a concept of Functional Airspace Blocks (FAB) is developed by Eurocontrol. A FAB is defined in the Single European Sky legislative package as "an airspace block based on operational requirements and established regardless European airspace is divided into nine Functional Airspace Blocks: NEFAB (North European FAB), Denmark-Sweden FAB, BALTIC FAB, FABEC (FAB Europe Central), FABCE (FAB Central Europe), DANUBE FAB, UK-IRELAND FAB and South West FAB.

of state boundaries, where the air navigation services and related functions is performance-driven and optimised through enhanced cooperation among ANSPs or an integrated provider" [START_REF]Single european sky legislative package[END_REF]. Nine FAB initiatives have been taken (Figure 1); two of these have already been implemented, namely the UK-Ireland and Denmark-Sweden FABs [START_REF]Eurocontrol Functional Airspace Block[END_REF].

Establishment of FABs will have consequences for NMOC operations. One of the key challenges will be to implement the centralized ATFM strategies in the individual FAB and develop the basis for the cooperation among the FABs [START_REF] Angenendt | DFS on Central Flow Managament Unit[END_REF]. Delegation of ATFM to individual FABs can lead to demand-capacity imbalances as each FAB may try to optimize traffic flow in its own airspace without considering other FAB requirements.

In this paper, we propose a distributed ATFM model which form the basis for interaction among FABs to implement ATFM strategies. The goal of distributed ATFM model is generate strategic 4D trajectories which minimize flight interaction in both three-dimensional space and in the time domain by extending a trajectory planning method developed by co-authors in [START_REF] Chaimatanan | Large Scale 4D Trajectory Planning[END_REF] and [START_REF] Islami | Large Scale 4D Trajectory Planning[END_REF].

The paper is organized as follows: Section II describes the concept of the proposed distributed ATFM model. Section III presents the proposed model and methodology in a mathematical framework. Then, a method to compute interaction between trajectories is presented in Section IV. Resolution algorithm to the problem is explained in Section V. Finally, numerical results are presented and discussed in Section VI.

II. CONCEPT DESCRIPTION

Centralized ATFM though offers a fair distribution and demand-capacity balance at a global level but stake-holders can only provide inputs and are not part of decision making. It may also lead to large number of pairwise reversals, i.e., the resulting order of flight arrivals can be quite different as compared to the original published flight schedules [START_REF] Bertsimas | Fairness in air traffic flow management[END_REF]. Because of this deviation from the original flight ordering, it may becomes difficult to implement such a solution locally.

In a distributed ATFM the decision-making responsibilities are shared between a number of airspace users (airlines, ANSPs, Military, Airports). Some examples of distributed ATFM are Ration by Schedule (RBS) [START_REF] Barnhart | Demand and capacity management in air transportation[END_REF] and Ground Delay Program (GDP) [START_REF] Charles | Stochastic optimization models for ground delay program planning with equity-efficiency tradeoffs[END_REF]. However the current state-of-the art is limited mostly to strategic planning, and the users' participation in planning reduces as the planning interval becomes smaller. With the gradual implementation of FAB concept in the European airspace, any distributed ATFM systems must also take into consideration multiple FAB interactions.

In a FAB scenario, the traffic flow management will be highly interdependent and will demand a significant cooperation with other FABs. One way to achieve this is by having distributed 4D trajectory planning amongst FABs which can ensures conflict-free trajectory for each aircraft. This will then translates into trajectory based operations where aircraft are required to fly a negotiated conflict-free trajectory through respective FABs.

A. Proposed Model

The main idea proposed in the paper is a distributed decision making model which can enable effective information sharing among the FABs for 4D trajectory planning which are conflict free. The goal of the proposed method is to separate a given set of aircraft trajectories in space and time domain by allocating an alternative fight plan (route and departure time) to each fight in a given FAB. Instead of trying to satisfy the capacity constraint, we focus on minimizing the interaction between trajectories. Where an interaction between trajectories occurs when two or more trajectories have an effect on each other; for instance, when trajectories occupy the same space at the same period of time. Therefore, contrary to the concept of conflict, the measurement of interaction does not only refer to the violation of minimum separation requirements. It also allows us to take into account other separation criteria such as minimum separation time between aircraft crossing at the same point.

The proposed concept is developed as follows: as illustrated in Figure 2, airspace users (airlines, airports, ANSPs etc.) submits relevant information (flight plan, slots, capacity) to NMOC, which then applies centralized traffic flow management strategies to match demand with capacity and other airspace constraints and generate revised flight plans. These flight plans are then used as an input to a FAB-Flight interaction Matrix. It is a 2D matrix which captures the flight interaction information between and within FABs. One dimension of the matrix is termed Controlling FABs and the other dimension is termed Intermediate FABs. A Controlling FAB is defined as a FAB where a given flight is originated or activated (in case of an enroute flight entering European airsapce), whereas an intermediate FAB is defined as a FAB through which a given flight traverse (over fly), terminates (lands) or exits the European airspace.

B. FAB-Flight Interaction Matrix

The FAB-Flight Interaction Matrix is developed as follows: as illustrated in Figure 4, for N FABs in a given airspace A, a 2D matrix of N rows and N columns is developed. The row vector of the matrix represents the number of flight interaction caused by a Controlling FAB C j in the Intermediate FABs I i for i = 1 to N. The column vector of the matrix represents the number of flight interaction caused by the Controlling FABs C j for j = 1 to N in an Intermediate FAB I i

I i = INT (C j , I A ) INT (C j , I B ) ... INT (C j , I i ) (1)
C j =     INT (C A , I i ) INT (C B , I i ) .... INT (C N , I i )     (2)
For example, as illustrated in Figure 4 Therefore, the total number of flight interactions U, in a given FAB i, can be given by summing the column vector:

U i = N ∑ j=1 INT (C j , I i ) (3)
The total flight interaction V in a given airspace A (which comprises of N FABs) can be given by

V = N ∑ j=1 U j (4)
This can be further normalized to find out the relative contribution of each FAB, in overall flight interactions, for a given airspace A.

U j rel = U j /V (5)
for j = 1 to N.

C. Traffic Flow Management Strategy

The distributed Traffic Flow Management strategy is developed as follows: first the Intermediate FAB with highest number of flight interaction is identified as a candidate FAB (equation 6).

FAB I i = MAX(U A ,U B ...,U N ) (6)
Then for the FAB I i , the Controlling FAB C j which generated highest number of flight interaction is identified (equation 7).

FAB C j = MAX(INT (C A , I i ), INT (C B , I i ), ..., INT (C N , I i )) (7)
The ATFM strategies (Space-Time separation) are then applied on randomly selected (fitness-proportional selection) flights in Controlling FAB C j . Flight interactions are recomputed given the revised flight plans and the FAB-Flight Interaction Matrix is updated. This process is repeated until the FAB-Flight Interaction Matrix is interaction free. Figure 5 illustrates the updated process, where the decision made by each FABs are evaluated by the optimization process. Then, the information of interaction based on FABs-flight interaction matrix is feed back to each FABs, which then make new decision and repeat the process until a solution that leads to minimum overall interaction is reached.

III. MATHEMATICAL MODELLING

This section sets the mathematical framework of the distributed air traffic flow management methodology we are proposing. First, a definition of interaction between trajectories is given. Then, the route / departure-time allocation techniques adapted for the distributed ATFM model is presented.

A. Interaction between trajectories

The concept of interaction between trajectories is introduced in [START_REF] Chaimatanan | Large Scale 4D Trajectory Planning[END_REF] and [START_REF] Islami | Large Scale 4D Trajectory Planning[END_REF]. It is a measurement that indicates when two or more trajectories occupy the same space at the same period of time. It is different from the conflict situation, which corresponds simply to a violation of the minimum separation (i.e. 5 NM horizontally and 1,000 ft vertically). Additional separation conditions, such as

! i,k = 2 d h < N h N h P i,k d h < N h
Fig. 6. Interactions, Φ i,k , at sampling point P i,k of trajectory i.

time separation, topology of trajectory intersection, distance between trajectories, etc. can also be taken into account in the concept of interaction. Consider a given set of N discretized 4D trajectories, where each trajectory i is a time sequence of 4D coordinates, P i,k (x i,k , y i,k , z i,k ,t i,k ), specifying that aircraft must arrive at a given point (x i,k , y i,k , z i,k ) at time t i,k , for k = 1, . . . , K i , and K i is the number of sampling points of trajectory i.

Consider a point k of trajectory i, interactions at point P i,k , denoted Φ i,k , may be defined as the total number of times that the protection volume around point P i,k is violated. Figure 6 illustrates interaction in the horizontal plane between N = 3 trajectories measured at point P i,k .

The interaction associated with trajectory i, denoted Φ i , is therefore defined to be:

Φ i = K i ∑ k=1 Φ i,k . (8)
Finally, the total interaction between trajectories, Φ tot , for a whole traffic situation is simply defined as:

Φ tot = N ∑ i=1 Φ i = N ∑ i=1 K i ∑ k=1 Φ i,k . (9)
One can observe that the measurement of the interaction between trajectories implicitly take into account the duration of conflict between trajectories. A practical methodology to compute the value of the interaction between trajectories in a large-scale context is presented in Section IV.

B. Route/Departure-time allocation

In order to separate the trajectories in 3D space and time domain, we rely on a route/departure-time allocation techniques introduced in [START_REF] Chaimatanan | Large Scale 4D Trajectory Planning[END_REF] and [START_REF] Islami | Large Scale 4D Trajectory Planning[END_REF]. The objective is to find alternative 4D trajectory for each flight, so as to minimize the total interactions between trajectories.

Given data. A problem instance is given by:

• A set of initial N discretized 4D trajectories with associated controlling FAB; • The discretization time step, ∆t;

• The number of allowed virtual waypoints, M;

∆ i :={-N i a .δ s , -(N i a -1).δ s , . . . , -δ s , 0, δ s , . . . , (N i d -1).δ s , N i d .δ s }. (10)
Alternative trajectory design. In this work, an alternative trajectory is constructed by placing a set of virtual waypoints, denoted

w i = {w m i |w m i = (w m ix , w m iy )} M m=1 , (11)
near the initial en-route segment and then by reconnecting the successive waypoints with straight-line segments as illustrated in Fig. 7. To limit the route length extension, the alternative en-route profile of flight i must satisfy:

L i (w i ) ≤ (1 + d i ), (12)
where L i (w i ) is the length of the alternative en-route profile determined by w i . Fig. 7 illustrated initial and alternative trajectories, constructed with M = 2 waypoints, where the location of each waypoint is constrained to be in a rectangular-shape possible location. Let W m ix be a set of all possible normalized longitudinal locations of the m th virtual waypoint on trajectory i. For each trajectory i, the normalized longitudinal component, w m ix , is set to lie in the interval:

W m ix := m 1 + M -b i , m 1 + M + b i , ( 13
)
where b i is a (user-defined) parameter that defines the range of possible normalized longitudinal component of the m th virtual waypoint on trajectory i. To obtain a regular trajectory, the normalized longitudinal component of two adjacent waypoints must not overlap, i.e.

m 1 + M + b i < m + 1 1 + M -b i (14)
and hence the user should choose b i so that

b i < 1 2(M + 1) . ( 15
)
Let W m iy be a set of all possible normalized lateral locations of the m th virtual waypoint on trajectory i. Similarly, the normalized lateral component, w m iy , is restricted to lie in the interval:

W m iy := [-a i , a i ], (16)
where 0 ≤ a i ≤ 1 is a (user-defined) model parameter that defines the range of possible normalize lateral location of the m th virtual waypoint on trajectory i, chosen a priori so as to satisfy (12). Let us set the compact vector notation: δ := (δ 1 , δ 2 , . . . , δ N ), and w := (w 1 , w 2 , . . . , w N ),

We shall denote by u i the components of u. It is a vector whose components are related to the modification of the i th trajectory, thereby our decision variable is:

u := (δ, w).
Finally, the interaction minimization problem can be formulated as a mixed-integer optimization problem, as follows:

δ i ∈ ∆ i , w m ix ∈ W m ix , w m iy ∈ W m iy , for all i = 1, . . . , N, m = 1, . . . , M, (P1)
where W m ix , and W m iy are defined by (13), and (16) respectively.

IV. INTERACTION DETECTION

In order to evaluate the objective function, at a candidate solution, (u), one needs to compute interaction between the N aircraft trajectories. To avoid the N(N-1) 2 timeconsuming pair-wise comparisons, which are prohibitive in our continental-scale application context, we propose a grid-based interaction detection scheme which is implemented in a so-called hash table as presented in [START_REF] Chaimatanan | Large Scale 4D Trajectory Planning[END_REF] and [START_REF] Islami | Large Scale 4D Trajectory Planning[END_REF].

First, the airspace is discretized using a four-dimensional grid (3D space + time), as illustrated in Figure 8. The size of each cell in the 4D grid is defined by the minimum separation requirement and the discretization time step, ∆t (see below). Then, for each given 4D coordinate P i,k (x i,k , y i,k , z i,k ,t i,k ) of each trajectory i, we identify which cell, says C i, j,k,t , of the 4D grid contains

P i,k (x i,k , y i,k , z i,k ,t i,k ).
Next, we consider each such cell C i, j,k,t and we successively check its surrounding cells (there are 3 3 = 27 such neighbouring cells, including cell C i, j,k,t itself). If one cell is occupied by an aircraft other than aircraft i itself, the horizontal distance (d h ) and the vertical distance (d v ) between the corresponding aircraft coordinates are measured. A violation of the protection volume is identified when both d h < N h and d v < N v .

In order not to underestimate interaction, and to avoid using small value to ∆t which leads to large number of trajectory samples and long computational time, we propose an inner-loop algorithm, detecting interaction between two sampling times, t and t + ∆t, by interpolating aircraft positions with a sufficiently small step size, t interp . The interpolation is performed only if no interaction is detected at time t. Then, one checks each pair of these interpolated points. The algorithm stops when an interaction is identified or when every pair of the interpolated points has been checked.

V. RESOLUTION ALGORITHMS

The 4D trajectory planning methodology for a distributed ATFM model, presented in this paper, relies on the interaction minimization problem introduced in Section III whose objective function values are obtained by simulation through the interaction detection scheme developed in Section IV. To solve the problem, a hybrid metaheuristic approach adapted to handle an air-traffic assignment problem at the continent scale is used. It relies on a classical simulated annealing (SA) algorithm and two different localsearch (LS) modules. The LS allows the system to intensify the search around a potential candidate solution while the SA allows the system to escape from a local trap and thereby ensuring the exploration of the solution space. The proposed hybrid algorithm combines the SA and the local search algorithm such that the LS is considered as an innerloop of the SA, which will be performed when a pre-defined condition is satisfied. The structure of the proposed hybrid algorithm of SA and LS methods is illustrated in Figure 9.

This hybrid SA-LS optimization algorithm has been applied to solve 4D trajectory planning at continent scale in [START_REF] Chaimatanan | Large Scale 4D Trajectory Planning[END_REF] and [START_REF] Islami | Large Scale 4D Trajectory Planning[END_REF]. In these works, the objective is, also, to minimize the total interaction between trajectories. However, the air traffic flow were managed based on a centralized decision making. In order to apply the hybrid SA-LS algorithm to our proposed distributed ATFM model, the LS search strategy and the neighborhood structure are modified according to our proposed ATFM strategy described in Subsection II-C (to be described below).

For our problem, the simulated annealing proceeds as follows. First, we evaluate the objective function at the current configuration (w, δ) C . It is denoted Φ C . Then a neighboring solution, (w, δ) N , is generated by a neighborhood function. Then, a new solution for this chosen flight is generated according to a pre-defined neighborhood structure. If the neighborhood solution improves the objective function value, then it is accepted. Otherwise, it is accepted with a probability e -∆Φ T , where ∆Φ = Φ N -Φ C is the difference of energy between current state C and new state N. When the maximum number of iterations, n T , at a given temperature is reached, the temperature is decreased according to the user-provided pre-defined schedule, and the process is repeated until the pre-defined final temperature, T f inal , is reached. More detail on simulated annealing can be found, for instance, in [START_REF] Dreo | Metaheuristics for hard optimization[END_REF].

Local search modules. The local search modules we use are heuristic methods that accepts a new solution only if it yields a decrease of the objective function. The process repeats until no further improvement can be found or until the maximum number of iterations n T LOC is reached. The two local-search modules correspond to the two following strategies:

• Intensification of the search on one Particular Trajectory (PT). Given a flight i, this state-exploitation step focuses on improving the current solution by applying a local change from the neighborhood structure only to flight i (only the decision variables (w i , δ i ) are affected).

• Intensification of the search on the Interacting Trajectories (IT). Given a flight i, this state-exploitation step applies a local change, from the neighborhood structure, to every flight that is both subjected to the same controlling FABs as flight i, and currently interacting with flight i.

Neighborhood structure. The hybrid algorithm we are proposing relies on a neighborhood structure to determine the next move. First, one has to determine which flight to be modified. In the framework of a distributed ATFM we are proposing, first one has to determine the controlling FABs which generates the highest proportion of interaction. Then, one chose one of the interacting flights in such FABs if Φ i ≥ τ • Φ avg , where τ is a user-defined parameter and Φ avg = Φ tot /N is the average value of interaction.

In order to generate a neighborhood solution for a given flight, i, from the current configuration (w i , δ i ) C , one has to determine whether to modify the location of waypoints or to modify the departure time in the next move. In general, searching for the solution in the time domain would be more preferable since it does not induce extra fuel consumption. However, empirical tests show that limiting the search to only that degree of freedom results in prohibitive computational time. Therefore, we introduce a user-defined parameter P w to control the probability to modify the location of the waypoints w i and such that the probability to modify rather the departure time is 1 -P w . For a given flight i, the neighborhood operator generates a new set of virtual waypoints or a new alternative departure time according to this probability P w .

Hybrid algorithm (SA and LS). Here is how the above-mentioned methods are combined. The methods are carried out according to pre-defined probabilities, which are proportional to the control temperature, T . The probability to carry out simulated annealing step, P SA , is:

P SA (T ) = P SA,min + (P SA,max -P SA,min ) • T 0 -T T 0 , (17)
where P SA,max and P SA,min are the maximum and minimum probabilities to perform the SA (pre-defined by the user).

The probability of running the LS module, P Loc , is given by:

P Loc (T ) = P Loc,min + (P Loc,max -P Loc,min ) • T 0 -T T 0 , (18)
where P Loc,max and P Loc,min are the maximum and minimum probabilities to perform the local search (defined analogously). And, finally the probability of carrying out both SA and the local search (successively), P SL , is:

P SL (T ) = 1 -(P SA (T ) + P Loc (T )). (19)
A key factor in tuning this hybrid algorithm is to reach a good trade off between exploration (diversification) and exploitation (intensification) of the solution space.

VI. NUMERICAL EXPERIMENTS

The proposed distributed ATFM model is implemented in Java. The overall methodology is tested with air traffic data involving flights over the European FABs, consisting of nine FABs listed in Table I. First, it is tested with a set of traffic consist of 4,000 flights over the European FABs on a UNIX platform with 1.7 GHz processor and 8 GB memory. The parameter values chosen to specify the optimization problem are given in Table II. The parameter values that specify the resolution algorithm are given in Table III. Then, it is tested with a full day en-route air traffic over the European FABs consisting of 26,122 flights on a UNIX platform with 2.4 GHz processor and 32 GB memory, using the same parameters as given in Table II and III, except that this time the number of iteration N I is set to 2,700. The initial and final total interaction between trajectories, the computation time, and the number of iterations performed to solve the distributed problems compared to the one solved based on centralized decision making methodology (proposed in [START_REF] Chaimatanan | Large Scale 4D Trajectory Planning[END_REF] and [START_REF] Islami | Large Scale 4D Trajectory Planning[END_REF]) are reported in Table IV. Evolution of the interaction between the 4,000 trajectories in each FAB and in the FAB-Flight Interaction Matrix during the optimization process using distributed and centralized models are presented in Figure 10, 11, 12, and 13 respectively. One can observe in Figure 10 and 12 that the FABs that has the highest level of initial interaction is FAB number 5 (FAB Europe Central). The controlling FABs which generates the highest interaction in FAB number 5 are FAB number 1,4, and 5. In the case of distributed model, the resolution algorithm tries to minimize the total interaction by modifying flight plans of flights associated to FAB number 1,4, and 5 before other FABs. While, in the centralized model case, the resolution algorithm modifies flight plans of flights that involves in Similarly, the evolution of the interaction between the full-day traffic, consisting of 26,122 trajectories, in each FAB and in the FAB-flight matrix during the optimization process using distributed and centralized models are presented in Figure 14, 15 , 16, and 17 respectively.

Again, one can observe that the FABs that in the case of distributed model, the resolution algorithm tries to minimize the total interaction by modifying flight plans of flights associated to the FABs that generate high interaction to the global system before other FABs. Although the trajectories can be separated only by modifying the horizontal flight profile and the departure time of each flight, the resolution algorithm finds an interaction-free solution in both Evolution of the FAB-flight interaction matrix for 26,122 trajectories during the optimization process using centralized ATFM model. Fig. 17. Evolution of the interaction in each FAB for 26,122 trajectories during the optimization process using centralized ATFM model. case 1 and case 2. In the case of small data set (case 1), the computation time to reach interaction-free solution are not significantly different when using distributed and centralized model. However, when the number of traffic increases (case 2), the distributed model converges to interaction-free solution faster than the centralized model. This is because the distributed ATFM strategy targets more on the FABs that generates high interaction. Nevertheless, both models yield interaction-free solution with computation time that is still viable for strategic planning and pre-tactical planning purpose.

VII. DISCUSSIONS

In this paper, we have presented a distributed air traffic flow management model, aiming to minimize total interaction between aircraft 4D trajectories at continent-scale. The objective was to develop a basis of information exchange and interaction among FABs for implementing distributed AFTM. The overall methodology is implemented and tested on a continent-size air traffic over the European FABs, and then compared with the one obtained using centralized ATFM model. The methodology based on distributed model converges to interaction-free solution faster than the centralized model. Furthermore, the color maps shows that the convergence rate of the distributed model is much better than centralized model and is viable for strategic planning as well as pre-tactical planning purpose. However, in practice, interfaces and different implementations, time delays and so forth need to be managed which present difficulties. These difficulties may not be compensated by these theoretical advantages.

Fig. 1 .

 1 Fig. 1.European airspace is divided into nine Functional Airspace Blocks: NEFAB (North European FAB), Denmark-Sweden FAB, BALTIC FAB, FABEC (FAB Europe Central), FABCE (FAB Central Europe), DANUBE FAB, UK-IRELAND FAB and South West FAB.

Fig. 2 .

 2 Fig. 2. Proposed Concept for information sharing among the FABs for interaction-free trajectory planning.

Fig. 3 .

 3 Fig. 3. An airspace divided into three FABs (FAB A, FAB B and FAB C) with four flight ( A, B, C and D) scenario.

Fig. 4 .

 4 Fig. 4. FAB-Flight Interaction Matrix

  , flight interactions due to flights controlled by FAB C A in the Intermediate FAB I A is given by row FAB C A and column FAB I A and denoted by INT (C A , I A ). Similarly, the number of flight interactions due to flights controlled by FAB C C in the same intermediate FAB i.e. FAB I c is given by row FAB C c and column FAB I c .

Fig. 5 .

 5 Fig. 5. FAB-Flight Interaction Matrix update process.

Fig. 7 .

 7 Fig. 7. Initial and alternative trajectories with rectangular-shape possible location of M = 2 virtual waypoints.

Fig. 8 .

 8 Fig.8. Four-dimension (3D space + time) grid for conflict detection, illustrated as a time series of 3D grids which is sampled with discretization time step ∆t = t nt n-1 . The size of each cell in the 3D grids is defined by the minimum separation requirements (N h and N v ).

Fig. 9 .

 9 Fig. 9. Structure of the proposed hybrid algorithm of simulated annealing and hill-climbing local search methods.

Fig. 10 .Fig. 11 .

 1011 Fig.10. Evolution of the FAB-flight interaction matrix for 4,000 trajectories during the optimization process using distributed ATFM model.

Fig. 12 .

 12 Fig. 12. Evolution of the FAB-flight interaction matrix for 4,000 trajectories during the optimization process using centralized ATFM model.

Fig. 13 .

 13 Fig.13. Evolution of the interaction in each FAB for 4,000 trajectories during the optimization process using centralized ATFM model.

Fig. 14 .

 14 Fig. 14. Evolution of the FAB-flight interaction matrix for 26,122 trajectories during the optimization process using distributed ATFM model.

Fig. 15 .

 15 Fig. 15. Evolution of the interaction in each FAB for 26,122 trajectories during the optimization process using distributed ATFM model.

Fig

  Fig. 16.Evolution of the FAB-flight interaction matrix for 26,122 trajectories during the optimization process using centralized ATFM model.

TABLE II CHOSEN

 II (USER-DEFINED) PARAMETER VALUES FOR THE OPTIMIZATION PROBLEM. Parameter Value Discretization time step, ∆t 20 seconds Discretization time step for possible departure-time shift, δs 20 seconds Maximum departure time shift, delta i a = δ i d := δ 120 minutes Maximum allowed route length extension coefficient, di 0.20 Maximum allowed flight level shifts, li,max := lmax 2 Maximum number of virtual waypoints, M 3

  Maximum probability to perform SA step, P SA,max 0.9 Minimum probability to perform local search step, PLoc,min 0.4 Maximum probability to perform local search step, PLoc,max 0.6 Number of iterations at each temperature step, NI 400 Number of iterations of the inner-loop local search step, nT LOC 5 Geometrical temperature reduction coefficient, β

	TABLE III
	EMPIRICALLY-SET (USER-DEFINED) PARAMETER VALUES OF THE
	RESOLUTION METHODOLOGY.
	Parameter	Value
	Minimum probability to perform SA step, P SA,min	0.8
		0.99
	Final temperature, T f	(1/500).T 0
	Inner-loop interpolation sampling time step, tinterp	5 seconds
	Probability to modify horizontal flight profile, Pw	1/3
	Probability to modify flight level, P l	1/3
	Threshold value, Φτ	0.5 Φavg

ACKNOWLEDGEMENT

Authors would like to thank Leilla Zerrouki from Eurocontrol for useful discussions on FAB design that greatly improved the manuscript.

AUTHORS BIOGRAPHY

Sameer Alam is a senior lecturer in aviation at the University of New South Wales, Canberra, Australia and a visiting researcher at ENAC (French Civil Aviation University). He obtained his BS degree in maths from University of Rajasthan, India (1994), M. Daniel Delahaye is Professor and head of the optimization team of MAIA laboratory at ENAC (French Civil Aviation University). He obtained his engineer degree from the ENAC and master of science in signal processing from the national polytechnic institute of Toulouse in 1991. He obtained his PhD in automatic control from the Aeronautic and Space National School in 1995 and did a post-doc at the Department of Aeronautics and Astronautics at MIT in 1996. He get his tenure in applied mathematics in 2012. He conducts research on mathematical optimization for airspace design (sector design, SID-STAR design).

Eric Feron graduated from Ecole Polytechnique, France, Ecole Normale Suprieure, France, and Stanford University, USA with the BS, MS and PhD, respectively. Eric Feron has held tenured faculty appointments at MIT (1993MIT ( -2005)), Georgia Tech (2005-Current), ONERA (1993)(1994)(1995)(1996)(1997)(1998), and ENAC (09/2011-Current). His current research interests include airport and airspace dynamics,UAVs and control theory. Currently, he is a Professor in the School of Aerospace Engineering at Georgia Tech.