Céline Pierdet
email: c_gedeon2000@yahoo.fr

FOR A CONTROL OF THE URBANISATION OF EMBANKMENTS IN RIVER TOWNS OF SOUTHEAST ASIA

Keywords: Urbanisation, riverfront, investor, Phnom Penh, Mekong river

Embankments in river towns of Southeast Asia were traditionally occupied by populations who practiced fishing and grew markets gardening. Since a few years we can observe a profound change of theses activities and of the urbanization. Traditional urbanisation is replaced with urban projects of great size, even with a continuous urbanisation which transforms them in waterfronts. So we have to try to understand why water seems to become today a vector of modernity and to create a recomposition of urban centrality round waterfronts. Some new urban projects are building in Phnom Penh South with important foreign capitals. Price of the soil is increasing on the embankments while most poor people are evicted from this attractive sector of the Cambodian capital.

So is this process common to others river towns in Southeast Asia such as Bangkok or Shanghai? For a long time were embankments in a neglected state because they were often insalubrious. Are they now destined to become common urban showcases for international investments?

urbain. L'urbanisation traditionnelle fait place à des projets de grande ampleur, voire à une urbanisation continue des berges, les transformant en waterfront. Dans quelle mesure l'eau ne constitue-t-elle pas aujourd'hui un vecteur de modernité et n'engendre-t-elle pas une recomposition de la centralité urbaine autour des fronts urbains fluviaux ?

De nouveaux projets d'urbanisme voient le jour à Phnom Penh au sud de la ville, caractérisés par un changement d'échelle par rapport à l'urbanisme traditionnel, et par des investissements de capitaux étrangers importants. Le prix du sol ne cesse d'augmenter sur ces terrains remblayés, provoquant l'éviction des populations les plus pauvres dans un secteur de la capitale redevenu très attractif.

Ce processus d'occupation des berges par un dense tissu de buildings s'est déjà produit à Bangkok ou à Shanghai. Est-ce dans cette logique que la capitale cambodgienne s'inscrit à son tour ? Les berges des métropoles asiatiques, longtemps délaissées car souvent jugées insalubres, sont-elles aujourd'hui destinées à devenir de banales vitrines pour les investisseurs internationaux ?

Mots-clefs :

Urbanisation, berge, investisseur, Phnom Penh, Mékong

Introduction:

For some years only, we are surprised by the speed with which the riversides of Phnom Penh change. This Southeast Asian capital of only 1,3 million inhabitants (2004), located in the confluence of the Mekong river, the Sap river and the Bassac river, began to highlight its riversides very recently. Until now, the embankments located in the centre of the capital were considered by the Municipality and the inhabitants to be insalubrious areas. These embankments were not valorised. Their lands fertilized every year by the flood of the Mekong river were left to the fishermen, to the populations which put them in cultures during the decrease of the water level.

However, this confluence area, named Chaktomuk, has a great symbolic dimension, as we have already could explain it 1 , due to the importance of water in Cambodian culture. This symbolic dimension besides constituted a driving force in 2001, when the Municipality interested in the peninsula. But since public equipments were accomplished on riversides, the prices of the soil increased and plans were changed in part.

It is this process which it is now a question of beginning to explain in this paper. Some new urban projects are building at the mouth of the Bassac river, with important foreign capitals. Price of the soil is increasing on the embankments while most poor people are evicted from this new attractive sector of the Cambodian capital. This process began in the most part of other Asian fluvial metropolitan areas for a long time. Traditional urbanisation is replaced with urban projects of great size, even with a continuous urbanisation which transforms them in waterfronts. So we have now to try to understand why water seems to become today a vector of modernity in Phnom Penh and how it creates a recomposition of urban centrality round riversides too, from other examples taken in Asia. Are the Phnom Penh's riversides now destined to become also common urban showcases for international investments? At first, it will be a question of describing fast the urban plans which affect riversides since a dozen years now, but with a much quicker rhythm lately. Then, we shall explain why urbanisation in a vast alluvial lowland, as that of Mekong river, requires the realization of polders behind the embankment, and in what the permanence of a city depends at first on the resistance of tiny dikes. Finally, it will be a question of analysing the urbanization process of the riversides in the Asian cities and 1 C. Pierdet, "La symbolique de l'eau dans la culture cambodgienne -Fête des eaux et projets urbains à Phnom Penh", in Géographie et Cultures, Paris : CNRS-L 'Harmattan, n° 56, juin 2006, pp. 5 -22. showing in what the recent changes put in danger the stability of cities and hydrosystems.

Urban plans at Chaktomuk area since 1993

After the leaving of the Vietnamese troops from Cambodia in 1989 and the lifting of embargo, Phnom Penh Municipality attributed a new importance to riverfronts and to the confluence area, in a context of national reconciliation. We can identify three stages in this process. This part of the city had never been developed. After the end of the flooding, the Municipality began to build a dike along the Bassac river. The objective of the Municipality was to use the silt accumulated at the mouth: so it built a polder on which to extend the city. After that, the Municipality, with the financial help of private companies, planned a public garden and a walk along the river, and roads.

The dike was extended from the Cambodiana Hotel to the Monivong bridge.

With this plan, the Municipality contributed to protect the city from the floods and to improve the life's conditions of the inhabitants. In a second time, it allowed to reorganize the land ownership along the river.

The integration of the Changvar peninsula in the city

In year 2000, the embankments of the peninsula were still occupied with agrarian populations who were farmers or fishermen, or with restaurants and karaoke.

During the floods, the inhabitants of the peninsula suffered serious damages because the back of the embankment was flooded. Moreover since the beginning of 1970s, the peninsula progressed to the south, up to the royal Palace, because of an accumulation of silts. This process increased too with the great flooding in year 2000.

A third plan, which was financed by the Municipality and by donations of private companies, consisted in protecting the inhabitants of the peninsula from the floods and in attracting the inhabitants of the city and the tourists2 . For symbolic reasons too, they tried to stop the progression of the peninsula in front of the royal Palace. Dikes and quays were built to protect the inhabitants from the floods. And public gardens were planned around a conference centre. They are now very seen frequently by the Cambodian families, especially in the evening.

So in this central part of the Cambodian capital, these plans allowed the inhabitants to appropriate the riverfronts and the river which were traditionally occupied by agrarian populations. And in the same time, the Municipality valorised this central place with public equipments and protected the city against the floods.

The problem of the urbanisation in the alluvial lowlands and in the deltas

The question is to understand how a lot of cities have been developed along so powerful fluvial organisms, most of the time from a tiny bank, when the natural pressures are maximal.

Numerous natural pressures

The cities which have been built in the alluvial lowlands are among the most threatened by the floods. And in the deltas, they have to cope with the risk of flooding or of change of the course of river. So it is a necessity to control this water which is omnipresent in the deltas, in particular in the deltas which are situated in Asia of Monsoons, because of the high densities of population and the dangerous rivers.

According to the geographer François Durand-Dastès, there is an "interaction between the high densities and the rice cultivation. It is an activity with a high productivity which allows a numerous population to live (…). Inversely it mobilizes an important hand of work."3 So the urbanization, in spite of the natural pressures, is a necessity to control and valorise the land. For example, in the delta of the Chao Phraya river, where high densities of population allowed the development of rice cultivation and of horticulture -and vice versa -, the city of Bangkok is protected against the floods only with a system of dikes and pumping stations. This metropolitan area is situated at an altitude of only 2 meters in 25 km of the sea.

The extension of a city from a bank

The embankments result from a natural process: when there is an exceptional flood, the rising of the river sets silts on the banks. The embankments, which surround the alluvial lowland, rise up flood after flood. In numerous alluvial lowlands, two embankments surround the minor bed of the river. They are rarely submerged, but on the other hand, the low lands located between the embankments and behind them are flooded each year during the period of high waters. So in the regions where they were no protections against the floods, the men fitted to this pressure and they built villages and roads on the embankments4 .

That's why the extension of a city from a bank, behind the embankment, requires to fill up the low lands and to protect the city from floods with a system of dikes. For example, in Europe, in the Netherlands, the abbeys developed polderisation process, with building of dikes and filling up the lowlands, from the Middle Ages, what favoured the urbanization. But in the majority of the countries, the polderisation process began only after the industrial revolution. So the geographer Jacques Bethemont explains that it allowed the transformation of "towns on water into towns on ground"5 .

On the other hand, from the natural embankments which surround the alluvial lowland and which protect several cities from floods in the same time, a second stage consisted in building quays in each city when the natural embankments were not enough effective. The quays protect the cities against the floods, make easy the road traffic and separate the urban activities from the fluvial activities.

We had to insist now on the fact that the natural embankments and the artificial dikes are tiny and long buildings which are very brittle.

The brittleness of the dikes and the embankments

In the deltas, the towns have the advantage of the proximity of the river and of the sea. Nevertheless the urban networks are not much developed. Most of the time, there is only one city which is located at the head of the delta because this site makes easier passage from a bank to the other one. In spite of the evolution of techniques, cities remain very vulnerable in these quasi amphibious areas. When the dikes break, the floods can be very impressive.

Recently, cases of breaks of dikes increased in the world everywhere. The case of News-Orleans, located in the delta of Mississippi, was very mediatized since the passage of cyclone Katrina on August 29th, 2005. This city settled in 1718 on a bank of this delta, under the level of waters of the lake Pontchartrain. Experts had established for a long time that resistance and height of dikes were deficient in case of cyclone of category 5. Besides, these dikes were very badly maintained. But their comments were not taken into account and the embankments of the lake Pontchartrain broke on several hundred of metres. The city was flooded under several metres of water. Other cases of breaks of dikes occurred for instance in April 2006 in the delta of Danube, following strong floods.

So it is really important to understand to what extent these tiny buildings are

brittle and require a constant maintenance, while urban pressure augments continuously, due to the acceleration of urban growth and attraction of the waterfronts.

The speeded up urbanisation process of the riversides in the Asian cities

The recent evolution registered on the riversides of Phnom Penh reflects the evolution registered in other Asian metropolitan areas for about twenty years. This ten-storeyed building exceeds broadly the height of buildings until then in force on the riversides of Phnom Penh, as well as the Monument of Independence. In the same way, in 2001, in an interview with the review Cambodge Nouveau, the manager of the Bureau des Affaires urbaines (BAU) explained that the Municipality planned to make dig the Pich island, or even to suppress it. Since 1980s, this island is cultivated as soon as the level of the river decreases. But the inhabitants of the island were evicted in 2005. The island was sold by the Municipality to a subsidiary company of Canadia Bank for the realization of a gated community, in exchange for the building of bridges on the Bassac river to ease congestion in the traffic of the bridge Monivong.

The recent evolution of the

So we can say that as soon as public equipments are planned or accomplished, the private capitals concentrate on the zone. As long as riversides were considered to be insalubrious areas, they were left in activities with weak added value.

With the increase in prices of the soil, plans are often changed in a speculative purpose, what produces a reconstruction of the urban centrality to the advantage of the Chaktomuk area.

The riversides of the Asian metropolitan areas: traditional urban showcases for the foreign investors, today in full mutation

In the most part of the fluvial metropolitan areas of Asia, the riversides support today high densities of buildings. But for instance, in Bangkok, the left bank of the river is occupied since the beginning of the 20th century by monasteries, but also by hotels, banks, a post office, the customhouse, etc.6 These ancient buildings, among which some of them still exist, are located just by the waterside. Riversides are almost completely privatised and this process is ancient. Today even, riversides are occupied by great hotels as the hotel Sheraton, or central post office. There is not public area to go for a walk along the river, on the left bank, in the centre of the city.

In the same way in Phnom Penh, it was impossible to accomplish a promenade all in one block along the river, going of the walk of the railway station to the Monivong bridge, due to hotels and other private buildings which advance to on the riverside.

The Cambodiana hotel interrupts the continuity of the walk. We can take another example with the city of Shanghai, which is located on the Huangpu river, in 30 km of the confluence with the Yangzijiang, and in 80 km of the sea. According to Pierre Gentelle7 , from 1930s, the towpath which edged Huangpu river was transformed into a walk, today named The Bund, along which being constructed about fifty buildings: seats of foreign banks, of foreign societies, etc. Therefore there also, the attraction for riversides is ancient. In April 1990, a new zone of business, named Pudong, was planned on the right bank of the Huangpu river. In The cities of Bangkok or Shanghai are opened to the foreigners for a long time.

In these cities-harbours, the riversides constituted the front-door of the city. They supported with it the first picture which presented itself to the visitors. And considering difficulties of urbanisation, the main commercial, administrative buildings, etc. were established on these riversides, or even industrial establishments. Today, in a very speculative context with a strong land pressure, the waterfronts know a very strong mutation with the establishment of numerous skyscrapers to increase the profitability.

It is now a question of re-qualifying these built-up areas buildings of which are now sometimes ancient. But Claude Chaline 8 , who broadly studied the evolution of waterfronts in cities-harbours, particularly in the United States, showed well how new plans tended to make the waterfronts into urban showcases, plated on agglomeration.

The plans are often fragmented, devoid of general vision, and left on private initiatives.

Most of the time, these plans are in dominant feature service sector, commercial, with some cultural activities. But they integrate badly with the whole agglomeration and with its functioning. Due to the presence of often futuristic architectural forms, as in Shanghai, which are located by the waterside, these waterfronts also become very attractive places, which it became necessary to visit. The geographer Chaline so a lot dwelled on the necessity to support pedestrian ways along riversides and to make them accessible to the inhabitants.

Urbanisation and hydrosystem in Southeast Asia: the increase of imbalances

At last, we had to understand that this increase of density of the soil on riversides or on these filled in areas has serious consequences to their stability. The development of metropolitan areas along these major rivers between others required the building of quays and dikes to protect them from dangerous floods. In the same way their development required the building of dams upstream to supply the inhabitants with electricity. But a river is a hydrosystem, what means that when we act on its course by constructing a dam, by accomplishing quays, etc it has consequences downstream. It changes above all the sedimentary balance sheet of the river, the processes of erosion and accumulation which are in work on embankments, but its flow too. 8 C. Chaline (dir.), Ces ports qui créèrent des villes. Paris : éd. L 'Harmattan, 1994, pp. 45-46. For example on the eastern embankment of the Changvar peninsula, the owners of the houses which are built at the edge of the Mekong river must build walls to protect themselves from floods, but also from the erosion which affects this bank (cf. photo 4). The articles of Cambodian newspapers relate more and more often that rice fields were devastated by a violent flood, or that riversides collapsed, drawing away in their fall some houses. Accumulations of silts bother the circulation or the fishing. All these problems are not new, but they have tendency to worsen. And especially, they move regularly.

It is necessary to understand that the embankments of a river are part of its hydrosystem. A too intensive urbanization of riversides risks worsen certain

1. 1

 1 Photo 1: Quay Sisowath, June 1992 (ECPAD -M. Riehl)

Photo 2 :

 2 Photo 2: The Changvar peninsula, June 2001 (Phnom Penh Municipality)

  Photo 3: The mouth of the Bassac river (C. Pierdet, September 2003)

  15 years, this quarter which was before occupied partly by factories, was absolutely transformed. The massive influx of State capitals, as well as those of the Municipality and of foreign investors allowed to hire development and incorporation speeded up of this quarter in the city. The quarter was linked to the old centre by the realization of two bridges, a tunnel and a line of subway. Today The Bund constitutes the central business district (CBD) of Shanghai. It concentrates a very big number of foreign banks, and the Stock Exchange of Shanghai. More than 200 skyscrapers in futuristic forms were constructed there, such as the Jinmao tower or the TV tower, which are famous all over the world.

Photo 4 :

 4 Photo 4: The eastern embankment of the Changvar peninsula, along Mekong river (C. Pierdet, September 2003)

The peninsula is linked up with Phnom Penh by a bridge which crosses the Sap river. The bridge was bombed in October 1972 and opened again to the circulation in 1993. Its reconstruction was financed by Japanese cooperation.

F. Durand-Dastès, Les systèmes d'utilisationde l'eau. CDU et SEDES, 1977, p. 67.

F. Durand-Dastès, op. cit., 1977, pp. 42-43.

J. Bethemont, "Villes et deltas", in La ville et le fleuve -Colloque tenu dans le cadre du 112 e congrès national des sociétés savantes(Lyon, 21-25 avril 1987). Paris : éd. duCTHS, 1989, pp. 330-331.

S.Van Beek, Bangkok, Then and Now. Bangkok: AB Publications, 2001, 131 p.

P. Gentelle, Chine et "Chinois" outre-mer : à l'orée du XXI e siècle. Paris :SEDES, 1999, p. 211.

imbalances, particularly for the populations who live on the embankments or in other cities situated downstream. It is brittle areas which ask for a frequent surveillance.

Conclusion:

These first elements of thinking allowed us to clear some stakes of the control of the urbanisation of riversides, particularly for the city of Phnom Penh. By comparison with plans accomplished in a metropolitan city as Shanghai, those of Phnom Penh are still of weak importance. However, the last plans of gated community in the mouth of Bassac seem already less well inserted into the whole metropolitan area. And the other plans have already had consequences to the evolution of all the riversides of the confluence area: some of them move back more and more, while silts continue settling elsewhere. It is necessary to encourage the Municipality to plan urban plans which are well inserted in the city and to preserve the public areas and the gardens which are located at the edge of the river and today approachable to all inhabitants. A general vision of the urban development must be favoured, and consequences of plans to the evolution of riversides must be anticipated and controlled. n° 56, juin 2006, pp. 5-22. Review Cambodge Nouveau, 2001. Van Beek S., Bangkok Then and Now. Bangkok : AB Publications, 2001, 131 p.

Bibliography


Emmanuel Pagnacco

E Souza De Cursi

R Sampaio2

Subspace inverse power method and Polynomial Chaos representation for the modal frequency responses of random mechanical systems

Keywords: Dynamic of structure, Frequency Response Function, Modal analysis, Propagation of uncertainty, Polynomial chaos, Clustered Eigenvalues, Rigid body modes

come

Subspace inverse power method and Polynomial Chaos representation for the modal frequency responses of random mechanical systems

Introduction

This study concerns the propagation of uncertainties for structural dynamic problems. Over the past decades, methods for solving stochastic eigenproblems or dynamics problems have been predominantly based on either Monte Carlo Simulation (MCS) or perturbation methods such as Neumann or improved perturbations methods [START_REF] Collins | The eigenvalue problem for structural systems with statistical properties[END_REF][START_REF] Hart | The treatment of randomness in finite element modelling[END_REF][START_REF] Shinozuka | Random eigenvalue problems in structural analysis[END_REF][START_REF] Dessombz | Analysis of stochastic structures: perturbation method and projection on homogeneous chaos[END_REF][START_REF] Szekely | Computational procedure for a fast calculation of eigenvectors and eigenvalues of structures with random properties[END_REF][START_REF] Schuëller | Uncertain linear systems in dynamics: Retrospective and recent developments by stochastic approaches[END_REF]. Implementation of MCS procedures is in principle straightforward while perturbation methods provide estimates of firsts moments, based on sensitivity analysis using Taylor series expansion of the quantity of interest. But these methods have several drawbacks. MCS is expensive in computing resources for large or complex problems or problems relying on several random variables, and perturbation methods have a limited radius of convergence for the involved series. Hence, efforts are constantly made to explore the suitability of spectral methods such as the Stochastic Reduced Basis Method (SRBM) [START_REF] Nair | New developments in computational stochastic mechanics, part ii applications[END_REF][START_REF] Nair | Stochastic reduced basis methods[END_REF] and the Polynomial Chaos (PC) representations [START_REF] Dessombz | Analysis of stochastic structures: perturbation method and projection on homogeneous chaos[END_REF][START_REF] Ghosh | A new algorithm for solving the random eigenvalue problem using polynomial chaos expansion[END_REF][START_REF] Ghanem | Efficient characterization of the random eigenvalue problem in a polynomial chaos decomposition[END_REF][START_REF] Panayirci | On the Capabilities of the Polynomial Chaos Expansion Method within SFE Analysis-An Overview[END_REF]. Both of them pertain to a non-statistical approach to represent randomness. SRBM is based on the subspace spanned by the considered application. In [START_REF] Nair | Stochastic reduced basis methods[END_REF], the practical basis suggested comes from a preconditioned stochastic Krylov subspace, performed from the nominal problem. Numerical studies for the first two moments on frequency response analysis of stochastic structural systems is addressed in this work, showing a degradation in results for large variations of the random variables. On the other hand, PC representation is a spectral representation of random processes which can be used in stochastic equations to represent unknown stochastic quantities [START_REF] Wiener | The homogeneous chaos[END_REF][START_REF] Ghanem | Stochastic Finite Elements A Spectral Approach[END_REF][START_REF] Maître | Spectral Methods for Uncertainty Quantification[END_REF][START_REF] Panayirci | On the Capabilities of the Polynomial Chaos Expansion Method within SFE Analysis-An Overview[END_REF]. Since it uses a Hilbertian framework, the PC expansion is able to represent any second-order random quantity with any probability law.

In this work, efforts are pursued to explore the suitability of spectral methods. More specifically, stochastic solutions in the frequency domain are sought for linear problems when random properties are considered in the mechanical model in the inertial and elastic properties. Evaluating these random results are generally referred to as Stochastic Finite Element Methods (SFEM). Linear dynamic response is e.g. more efficiently dealt with in the modal space. To deal with the modal space one has first to solve an eigenvalue problem.

This problem can either be a differential eigenvalue problem or a matrix eigenvalue problem, depending on whether a continuous model or a discrete model is used to describe the given system. For stochastic problems, probabilistic characterization of the eigensolutions of random matrix and of differential operators turn out to be of interest. SFEM with PC representation are well developed for linear algebraic systems, but are less developed for the random eigenvalue and for the random response of dynamic problems. The random eigenanalysis problem appears to be non linear and challenging conditions arise when part of the spectrum is clustered, that is, when there are repeated or closely spaced eigenvalues. In Ref. [START_REF] Dessombz | Analysis of stochastic structures: perturbation method and projection on homogeneous chaos[END_REF], the PC eigenvalue problem is approximated through a projection onto the deterministic normal mode basis, both for the normal mode equilibrium equation and for the normalization equation. In Ref. [START_REF] Ghosh | A new algorithm for solving the random eigenvalue problem using polynomial chaos expansion[END_REF][START_REF] Ghanem | Efficient characterization of the random eigenvalue problem in a polynomial chaos decomposition[END_REF], the eigenvalue problem with distinct roots is rewritten as a set of non-linear equations with a L2 normalization of the normal modes.

Next, specific strategies are developed in Ref. [START_REF] Ghosh | Stochastic convergence acceleration through basis enrichment of polynomial chaos expansions[END_REF][START_REF] Ghosh | An invariant subspace-based approach to the random eigenvalue problem of systems with clustered spectrum[END_REF] to deal with repeated roots when such an approach is carried out. However, while the linear stochastic dynamic problem is carried out without approximation in this way -apart from the practical necessity of polynomial truncature-, the involved computational effort appears to be much more important than the one involved for static random problems. Hence, to be more efficient numerically, the Ref. [START_REF] Verhoosel | Iterative solution of the random eigenvalue problem with application to spectral stochastic finite element systems[END_REF] proposes to transpose in the random domain the deterministic inverse power method based on the spectral shift strategy. It appears to be a robust method to determine the random eigenvector associated to a random eigenvalue. But normalization of normal modes is achieved only in the mean sense in this work. This does not lead to favorable conditions to address next modal response of dynamic problems.

In this work, stochastic discrete linear systems, or discretized continuous systems, having -or not-null or repeated eigenvalues are considered for obtaining modal based frequency responses. PC representation through a weak characterization of mechanical equations is focused and efforts are made to solve the resulting stochastic modal equations and frequency responses without transforming this problem as it would be if, for example, projections onto the deterministic normal modes basis are involved. To achieve a numerical efficiency, the Newton-Raphson strategy is not selected and the deterministic eigenanalysis method such as the deflated inverse power method and the subspace inverse power method are transposed to the stochastic framework. The modal problem is addressed here without approximations on the normalization condition for the normal modes, by requiring that the stochastic normalization and the stochastic orthogonalization conditions must be satisfied. Then, modal based frequency responses take advantage of the stochastic modal mass-orthonormalization which is ensured. Moreover, for numerical efficiency, it is also a main objective of this work to keep the computational costs similar to the one required for the -well known and used-stochastic static problem solved by the PC method. Indeed, one numerical drawback of the intrusive PC strategy is the use of expectations of multiple PC products which requires both lot of operations to be computed and lot of memory to be stored. This aspect is a part of the curse of dimensionality problem which arises for stochastic modal based frequency responses analysis and it must keep in mind, as we will see in the following of this paper that expectations of a five PC products could be required to address these problems (while only expectations of a three PC products are required to address stochastic static problems). This problem can occur both for ensuring orthonormalization of the modal base and for obtaining modal based frequency responses. Hence, we will pay special attention to propose alternative numerical strategies at both these steps in order to use only expectations of PC products which are already used for stochastic static problems. Adoption of these propositions ensure the efficiency of the complete stochastic modal based frequency response analysis of this work by having computational costs similar to the one required for the stochastic static problem solved by the PC method.

The paper is organized as follows: the Sect. 2 addresses the modal dynamic problem, considering the deterministic and the stochastic frameworks. The Sect. 3 recalls basic considerations of the PC representation and uncertainties propagation in the context of the dynamic responses, either from the direct or from the modal way. The Sect. 4 addresses the proposed numerical strategies for the inverse power method and its deflated and subspace variants, by transposing the deterministic strategies to the stochastic ones. Then, four application examples are chosen to demonstrate the applicability and effectiveness of the proposed methodologies in Sect. 5. Applications to eigenanalysis random problems as well as frequency modal random responses are carried out for three simple application examples and one realistic application example. It is then demonstrated that the proposed strategies are applicable and that the proposed method is able to deal with situations for which part of the spectrum is clustered. Finally, a summary is given in Sect. 6 to conclude this work.

Problem statement

Deterministic problem

Let us considers the discrete -or discretized-linear deterministic problem of computing the dynamic response of a system subjected to a forcing which is defined in the frequency domain (ω = 2πf ):

k -ω 2 m + jcω u (ω) = f (ω) (1)
where upright bold letters denote structural deterministic matrices. In this expression k, m, and c are ndimensional real symmetric structural matrices, being respectively the stiffness, mass, and damping matrices.

The frequency response functions are u (ω) when f (ω) is chosen as the n × n identity matrix i over a specified frequency band, ω min ≤ ω ≤ ω max .

Since the problem is linear with respect to the frequency response, it is interesting to introduce the modal representation associated to the above spatial representation. Restricting ourselves to proportional-damped systems, the modal representation is real and is defined as:

     k φ = m φ λ φ T m φ = µ (2)
with ξ = φ T c φ the modal damp matrix and where:

• φ is the normal-mode matrix formed with the eigenvectors as columns;

• λ is the diagonal matrix formed with the corresponding eigenvalues, that are the squares of the natural circular frequencies; and

• µ is the diagonal modal mass matrix that is imposed to be the identity matrix, i.e. µ = i, in this study.

Notice that a widely used convention consists in ordering the eigenvalues in ascending values, that is

λ 1 ≤ λ 2 ≤ • • • ≤ λ n ,
while the associated eigenvectors are:

φ = {φ 1 , φ 2 , • • • , φ n } .
Hence, when there is clustered eigenvalues, switching of the ordering of eigenvectors can occur for a small change in the physical parameters. For example, considering distinct natural frequencies of a simple plate with varying dimensions, ordering of a torsional and a flexural mode can be switched due to the modification of these frequencies. However, these two modes can be recognized experimentally, when compared before and after modifications, because they keep their shapes. Thus, a Modal Assurance Criterion (MAC) will indicate that they are nearly parallel, showing then that they correspond to the same physical modes1 .

Then, the solution u (ω) is such that:

u (ω) = φ x (ω) = n r=1 φ i x i (ω) (3)
for x (ω) which satisfies the uncoupled modal equations:

λ -ω 2 i + jξω x (ω) = φ T f (ω) .
However, for mechanical applications which involves many degrees of freedom, an approximate solution is generally sought practically by truncating the modal basis to a specified range where only the lowest eigenvalues separated from the rest of the spectrum are considered:

λ 1 ≤ • • • ≤ λ nr < λ nr+1 ≤ • • • λ n .
This leads us to rewritten Eq. 3 as:

u (ω) nr r=1 φ i x i (ω) = φ x (ω)
where φ becomes now a n × n r matrix, while modal matrices are adjusted in their sizes accordingly.

Stochastic problem

When the system has random parameters, its response becomes a stochastic process, U (ω, ), such that

K ( ) -ω 2 M ( ) + jC ( ) ω U (ω, ) = f (ω) (4)
for all ∈ Ω where (Ω, A, Prob) is the probability space associated with the underlying physical experiments and where upper case upright letters denotes random variables corresponding to deterministic ones which are denote by the same letter in lower case. The stochastic eigenproblem associated to this model is:

     K Φ = M Φ Λ Φ T M Φ = µ (5)
where the argument is dropped for brevity. At this point, one has to note that this modal problem is slightly different from the one of Ref. [START_REF] Verhoosel | Iterative solution of the random eigenvalue problem with application to spectral stochastic finite element systems[END_REF][START_REF] Ghanem | Efficient characterization of the random eigenvalue problem in a polynomial chaos decomposition[END_REF][START_REF] Ghosh | Stochastic convergence acceleration through basis enrichment of polynomial chaos expansions[END_REF][START_REF] Ghosh | An invariant subspace-based approach to the random eigenvalue problem of systems with clustered spectrum[END_REF] due to the requested mass-orthonormalization condition.

Then, U (ω) can be express as:

U (ω) = Φ X (ω)
for:

Λ -ω 2 i + jΞω X (ω) = Φ T f (ω)
where:

Λ = Φ T K Φ and Ξ = Φ T C Φ.
One can notice that, due to the possible switching, or even mixing, of eigenvectors in the deterministic problem, the stochastic eigenproblem can be difficult to solve when there is clustered eigenvalues [START_REF] Ghosh | An invariant subspace-based approach to the random eigenvalue problem of systems with clustered spectrum[END_REF].

3 Propagation of uncertainties using polynomial chaos representation Several methods can be applied to get a representation of U(ω) for given distributions of the random input properties [START_REF] Schuëller | Uncertain linear systems in dynamics: Retrospective and recent developments by stochastic approaches[END_REF][START_REF] Panayirci | On the Capabilities of the Polynomial Chaos Expansion Method within SFE Analysis-An Overview[END_REF]. From them, as an attractive alternative to the MCS, the PC representation consists in the expansion of the random process over a set of orthogonal polynomials. Wiener first introduced it in 1938 [START_REF] Wiener | The homogeneous chaos[END_REF] and used only Hermite polynomial basis to represent Gaussian random processes. More recently, works on generalized Polynomial Chaos propose a list of polynomial basis corresponding to non-Gaussian distributions [START_REF] Xiu | The Wiener-Askey polynomial chaos for stochastic differential equations[END_REF]. In effect, one can use any distribution as kernel and orthogonalize a polynomial basis of L 2 with respect to the probability measure generated by the chosen distribution. [START_REF] Ernst | On the convergence of generalized polynomial chaos expansions[END_REF] demonstrated the L 2 -convergence for any L 2 functionals, that is for any second-order random process (i.e. with finite variance).

Direct dynamic frequency response

Considering a generic random vector or random matrix, G, the PC representation proposes to express it as a polynomial (truncation of a series) using a set of n g + 1 orthogonal polynomials, denoted ψ g , in variables

ζ i , i ∈ {1, , 2, . . . d}, that are collected in a d-dimensional vector ζ : G( ) = ng g=0 Ψ g (ζ( ))g g (6)
where n g + 1 = (d+h)! d! h! , h being the order of the representation. Notice that italic bold letters are used for vectors of the spectral representation. Generally, for a practical implementation, the order of expansion results from a truncation which has to be chosen accordingly a suitable criterion. For the structural matrices

K (ζ), M (ζ), or C (ζ)
and for responses U(ω, ζ) the adopted PC representation are:

K (ζ) = n k k=0 Ψ k (ζ) k k , M (ζ) = nm m=0 Ψ m (ζ) m m , C (ζ) = nc c=0 Ψ c (ζ) c c ,
and:

U(ω, ζ) = nu u=0 Ψ u u u (ω) .
Coming back to the generic random vector or matrix G, the n g vectors or matrices of deterministic coefficients g g , now used to describe G, can be evaluated in two ways: using an intrusive or a non-intrusive method. The non-intrusive method can be always used to find the coefficients of the stochastic variables when they are represented by a PC expansion. It uses the orthogonality of the polynomials with respect to the appropriated inner product to evaluate each vector of deterministic coefficients g g :

g g = < G (ζ) , Ψ g (ζ) > < Ψ g (ζ) , Ψ g (ζ) > (7)
where

< G (ζ) , Ψ g (ζ) >= ´Rd G (ζ) Ψ g (ζ) p ζ (ζ)
dζ represents the inner product introduced by the PDF

p ζ (ζ).
To evaluate this inner product, it is introduced a partition, ζ i , the integrations points, w the vector which collects the quadrature weights and n gp the number of points, the above integration can be done using a Gauss quadrature rule ( [START_REF] Press | Numerical Recipes in C: The Art of Scientific Computing[END_REF]):

< G (ζ) , Ψ g (ζ) >= n GP i=1 G (ζ i ) Ψ g (ζ i ) {w} i .
In this way, the residual associated to the difference between G( ) and its PC representation is orthogonal to its PC representation. This strategy to determine the representation can be referred to as a strong characterization [START_REF] Ghosh | An invariant subspace-based approach to the random eigenvalue problem of systems with clustered spectrum[END_REF] and can be used to find coefficients of input stochastic variables such as k k , m m or c c for the structural matrices as well as to model the output stochastic variables such as u u (ω) for responses

U(ω, ζ).
The intrusive method, or weak characterization, follows a variational approach. As a first step, the PC representation introduced for the stochastic variables are put into the governing equations. For instance, Eq.

(4) produces:

n k k=0 Ψ k k k -ω 2 nm m=0 Ψ m m m + jω nc c=0 Ψ c c c nu u=0 Ψ u u u (ω) = f (ω) . (8)
Next, the vectors coefficients are given by solving the system which is obtained by taking expectations of the equations produced by the projection of Eq. ( 8) onto the set of polynomials Ψ:

nu u=0 n k k=0 Ψ i Ψ k Ψ u k k -ω 2 nm m=0 Ψ i Ψ m Ψ u m m + jω nc c=0 Ψ i Ψ c Ψ u c c u u (ω) = Ψ i f (ω) (9)
for i = 0, . . . , n u . This leads to a deterministic matricial system of the form:

K -ω 2 M + jCω U (ω) = F (ω) (10)
where:

K =          n k k=0 Ψ 0 Ψ k Ψ 0 k k n k k=0 Ψ 0 Ψ k Ψ 1 k k • • • n k k=0 Ψ 1 Ψ k Ψ 0 k k n k k=0 Ψ 1 Ψ k Ψ 1 k k . . . . . .          , U =       u 0 u 1 . . .       and F =       Ψ 0 f Ψ 1 f . . .      
.

where calligraphic letters denote the structural PC representation matrices used for the analysis. The matrices M and C have similar form as the matrix K. Equation ( 10) is a linear system of deterministic equations which involves expectations of three PC products.

From a numerical point of view, the weak characterization provides a set of n × n g coupled algebraic equations, while the strong characterization determines g g from integrations, requiring at least a set of experiments for a quadrature rule. When the order of representation is adequate, the coefficients obtained from both characterizations should be identical when the same orthogonal basis is considered. However, this is not true for an inadequate truncation since the weak characterization couples the coefficients, while the strong characterization determines each coefficient independently.

Modal frequency response

The direct strategy proposed in the previous section by using the intrusive method is expensive numerically when considering a large frequency band since it involves solving a large system many times. As the structural matrices do not depend on the frequency, the modal response presented in subsection 2.2 is relevant.

Considering only a particular, single mode, Φ r , modal quantities can be represented in the chosen basis as:

Φ r (ζ) = n λ φ=0 Ψ φ (ζ) φ rφ , Λ r (ζ) = n λ λ=0 Ψ λ (ζ) λ rλ
and the Eq. ( 5) , taken in the weak sense, becomes:

                     n λ φ=0 n k k=0 Ψ i Ψ k Ψ φ k k φ rφ = n λ φ=0 n λ λ=0 nm m=0 Ψ i Ψ m Ψ φ Ψ λ m m φ rφ λ rλ n λ φ=0 nm m=0 n λ λ=0 Ψ i Ψ φ Ψ m Ψ λ φ T rφ m m φ rλ = Ψ i µ (11)
for i = 0, . . . , n λ which can be rewritten in a matricial form as:

     K Φ r = M Φ r Λ r Φ T r M Φ r = Υ (12)
when collecting:

Φ r =       φ r0 φ r1 . . .       , Λ r =       λ r0 λ r1 . . .       , Υ =       Ψ 0 µ Ψ 1 µ . . .       =       Ψ 0 i 0 . . .      
and:

M Φ r =           n λ φ=0 nm m=0 Ψ 0 Ψ k Ψ φ Ψ 0 mmφ rφ n λ φ=0 nm m=0 Ψ 0 Ψ k Ψ φ Ψ 1 mmφ rφ • • • n λ φ=0 nm m=0 Ψ 1 Ψ k Ψ φ Ψ 0 mmφ rφ n λ φ=0 nm m=0 Ψ 1 Ψ k Ψ φ Ψ 1 mmφ rφ . . . . . .          
. This is a non-linear system of deterministic equations which involves expectations of four PC products. The reference [START_REF] Ghanem | Efficient characterization of the random eigenvalue problem in a polynomial chaos decomposition[END_REF] proposes to solve it using a Newton-Raphson method starting from a statistically sampled initial value.

Then, synthesis of U (ω) leads to consider:

U (ω) = nr r=1 n λ φ=0 Ψ φ φ rφ X r (ω) (13)
for X r (ω) = nx x=0 Ψ x x rx (ω) such that:

Φ T r KΦ r -ω 2 Φ T r MΦ r + jΦ T r C r Φ r ω X r (ω) = Φ T r f (ω)
or:

Φ T r KΦ r -ω 2 + jΦ T r CΦ r ω X r (ω) = Φ T r f (ω)
since normal modes are mass normalized. It is, in the weak sense:

nx x=0   n λ φ=0 n k k=0 n λ λ=0 Ψ i Ψ φ Ψ k Ψ λ Ψ x φ T rφ (k k + jωc k ) φ rλ -ω 2 Ψ i Ψ x x rx (ω) = Ψ 2 i φ T ri f (ω) (14)
for i = 0, . . . , n x . This is a linear system of deterministic equations which involves expectations of five PC products. Then it is more efficient computationally to use PC representation of modal quantities. Using

Λ r = n λ λ=0 Ψ λ λ rλ and Ξ r = n ξ ξ=0
Ψ ξ ξ rξ , the above expression is rewritten:

Λ r -ω 2 + jΞ r ω X r (ω) = Φ T r f (ω)
to produce, in the weak sense:

nx x=0 n λ λ=0 Ψ i Ψ λ Ψ x λ rλ -ω 2 Ψ i Ψ x +jω n ξ ξ=0 Ψ i Ψ ξ Ψ x ξ rξ   x rx (ω) = Ψ 2 i φ T ri f (ω)
for i = 0, . . . , n x . It is a linear system of deterministic equations which involves expectations of three PC products. This leads to a deterministic matricial system of the form:

J r -ω 2 I r + jD r ω X r (ω) = R r (ω) (15)
where:

J r =          n λ λ=0 Ψ 0 Ψ λ Ψ 0 λ rλ n λ λ=0 Ψ 0 Ψ λ Ψ 1 λ rλ • • • n λ λ=0 Ψ 1 Ψ λ Ψ 0 λ rλ n λ λ=0 Ψ 1 Ψ λ Ψ 1 λ rλ . . . . . .          , X r =       x r0 x r1 . . .       and R r =       Ψ 2 0 φ T r0 f Ψ 2 1 φ T r1 f . . .      
, while matrices I r and D r have a similar form to the matrix J r . Thus, to evaluate U (ω) through Eq. ( 13) once the stochastic eigenproblem is solved, assuming it provides the PC representation of eigenvectors and eigenvalues, it is efficient to solve Eq. ( 15) to obtain the modal displacement matrix coefficients (that is a vector, if f (ω) also is). It corresponds to the equation of a single degree of freedom system, thanks to the normal modes basis that decouples mechanical degrees of freedoms. This is the intrusive strategy which is chosen in this work to solve the stochastic dynamic problem by the modal strategy. But it remains necessary to propose an efficient numerical strategy to solve the system [START_REF] Wiener | The homogeneous chaos[END_REF] which defines the normal modes base. This is the subject of the next section.

Numerical methods for the determination of normal modes bases

Considering the deterministic generalized normal modes basis problem, it exists numerical methods dedicated to its efficient determination [START_REF] Géradin | Mechanical vibrations: Theory and applications to structural dynamics[END_REF][START_REF] Arbenz | Lecture notes on solving large scale eigenvalue problems[END_REF]. Let us first describe the ones we propose to transpose in the random domain.

Deterministic methods

From a numerical point of view, Lanczos method and power method are adapted to find eigenvectors of large numerical problems. The inverse power method -or inverse iterations method-is the most basic method of computing an eigenvector φ 1 associated to the lowest eigenvalue |λ 1 | > 0 for the generalized problem [START_REF] Hart | The treatment of randomness in finite element modelling[END_REF] when considering a positive-definite symmetric and real stiffness matrix. The basic algorithm of this method is:

1. Choose an initial vector φ (0)

1 such that φ (0)T 1 m φ (0) 1 = 1, 2. For k = 1, 2, . . . do: (a) Evaluate f (k) 1 = m φ (k-1) 1 (b) Solve k φ (k) 1 = f (k) 1 (c) Normalize φ (k) 1 such that φ (k)T 1 m φ (k) 1 = 1 3. Let φ 1 ← φ (k) 1 and compute λ 1 = φ T 1 k φ 1

End

This algorithm repeats the step 2 until φ

(k) 1
converges to within some tolerance . A variant is obtained by including the step 3 in the loop of step 2, enabling to test the convergence on λ

(k) 1 , being λ (k) 1 -λ (k-1) 1 < λ (k) 1 . Noticed that even if φ (0)
1 is orthogonal to φ 1 , it is expected that φ 1 would be recovered numerically due to round-off errors.

Then, deflation can be carried out to construct the next lowest eigenvalue λ 2 and so on for the next ones.

In such a way, the eigenvalue λ r , for r > 1, is obtained by introducing an orthogonalization step before the normalization, as:

φ (k) r = φ (k) r - r-1 q=1 φ q φ T q mφ (k) r
This orthogonalization step is also necessary when considering semi-definite stiffness matrix where the first n 0 eigenvectors φ RBM correspond to the null space of k, being kφ RBM = 0. In such situation, a generalized inverse is applied within iterations. Then, for n r elastic modes, the adopted deflated inverse power method algorithm is:

1. Find the null space φ RBM of k such that φ T RBM m φ RBM = i 0 and let φ ← φ RBM 2. For each sought mode r = n 0 + 1, n 0 + 2, . . . , n r , do:

(a) Choose an initial vector φ (0) r such that

φ T RBM m φ (0) r = 0 and φ (0)T r m φ (0) r = 1, (b) Compute λ (0) r = φ (0)T r k φ (0) r (c) For k = 1, 2, . . . do: i. Evaluate f (k) r = m φ (k-1) r ii. Solve k φ (k) r = f (k) r iii. Orthonormalize φ (k) r such that φ T m φ (k) r = 0 and φ (k) r m = 1 iv. Compute λ (k) r = φ (k)T r k φ (k) r (d) Let φ ← φ, φ (k) r 3.
End where i 0 is the n 0 identity matrix. However, this algorithm does not ensure to find the first eigenvectors when repeated eigenvalues exist. The inverse simultaneous vector iteration -or subspace inverse power methodhas not this drawback. It consists in handling simultaneously a set of modes φ (k) , while preserving their orthonormality within iterations. Its basic algorithm is:

1. Find the null space φ RBM of k such that φ T RBM m φ RBM = i 0 2. Choose an orthonormalized initial matrix φ (0) such that:

φ T RBM m φ (0) = 0 and φ (0)T m φ (0) = i r0 3. Compute λ (0) = φ (0)T k φ (0) 4. For k = 1, 2, . . . do: (a) Evaluate f (k) = m φ (k-1) (b) Solve k φ (k) = f (k)
(c) Orthonormalize φ (k) such that:

φ T 0 m φ (0) = 0 and φ (k)T m φ (k) = i r0 (d) Compute λ (k) = φ (k)T k φ (k) 5. Let φ ← φ RBM , φ (k)

End

where i r0 is the n r -n 0 identity matrix and φ (0) is the n × (n r -n 0 ) matrix composed of starting guesses vectors. If it is not deficient in the directions of sought eigenvectors, convergence of this method is ensured for real symmetric semi-definite positive matrix k.

4.2 Determination of the stochastic normal modes basis using the deflated inverse power method and polynomial chaos representation

The deflated inverse power method is a suitable alternative to the Newton-Raphson method to solve the system of Eqs. [START_REF] Panayirci | On the Capabilities of the Polynomial Chaos Expansion Method within SFE Analysis-An Overview[END_REF] in order to build the stochastic normal mode basis.

The proposed algorithm for the deflated inverse power method is the following:

1. Choose Φ RBM from null (k) to represent Φ RBM 2. Normalize Φ RBM such that Φ T RBM M Φ RBM = i 0 and let Φ ← Φ RBM 3.
For each sought mode r, do:

(a) From the deterministic vector φ r , define the initial vector Φ

r such that Φ T MΦ (k) r = 0 and let Φ (0) r ← Φ (0)1 r (c) Evaluate Λ (0) r from a modal projection of K over Φ (0) r (d) For k = 1, 2, . . . do: i. Evaluate F (k) r = M Φ (k-1) r ii. Solve K Φ (k) r = F (k) r iii. Find Φ (k)⊥ r that orthogonalizes Φ (k) r to Φ q for q ∈ {1, . . . , r -1} and let Φ (k) r ← Φ (k)⊥ r iv. Evaluate Φ (k)1 r
that normalizes modal mass to unity and let

Φ (k) r ← Φ (k)1 r v. Evaluate Λ (k) r from a modal projection of K over Φ (k) r (e) Let Φ ← Φ, Φ (k) r

End

In this algorithm, a simple and efficient starting guess for

Φ (0)
r is produced from the normalized normal mode φ r of the deterministic nominal mechanical system, in which case Φ (0) r0 = φ r while Φ (0) rφ = 0 for φ > 0. Hence, for an efficient numerical scheme, it is recommended to systemically carried out the deterministic subspace inverse power method of the previous subsection 4.1 in order to produce efficiently all the starting guesses Φ (0) r for r ∈ {n r+1 , . . . , n}. This choice ensure to find a meaningful set of first orthogonal eigenvectors, even when repeated eigenvalues exist for the nominal system. From our numerical experiment, it is found that this choice for the starting guess enables the algorithm to deal with deterministic nominal systems having repeated eigenfrequencies.

However, it has to be noted that even when the deterministic eigenvector is effectively mass normalized, its stochastic counterpart may not be mass normalized. Then, the normalization sub-step (b) of this algorithm holds in any cases. Modal projections, orthogonalization, and normalization steps are detailed in the following subsections. Notice that choosing the suggested starting guess for Φ (0) r enables the algorithm to deal with null eigenvalues, in which case the sub-step (d) is bypassed.

Normalization step

The normalization step consists in searching Φ (k)1 r which is defined by:

Φ (k)1T r M Φ (k)1 r = 1.
It is achieved for: where:

Φ (k)1 r = 1 Φ (k) r M Φ (k)
Φ (k) r M = Φ (k)T r M Φ (k) r 1 /2 while Φ (k)1 r (ζ) = n λ φ=0 Ψ φ (ζ) φ (k)1 rφ with φ (k)1 rφ = < Φ (k)1 r , Ψ φ (ζ) > < Ψ 2 φ (ζ) > .
In practice, it can be performed from a quadrature rule:

φ (k)1 rφ = 1 < Ψ 2 φ (ζ) > ngp i=1 Φ (k)1 r (ζ i ) Ψ φ (ζ i ) w i
where:

Φ (k)1 r (ζ i ) = Φ (k)T r (ζ i ) M (ζ i ) Φ (k) r (ζ i ) -1 /2 Φ (k) r (ζ i ) .
This is an enhanced formulation from a numerical point of view since:

• it does not involved the non-linear solution of the equations as the ones obtained by the weak characterization, cf. Eq. ( 11);

• it does not requires expectations of a four PC products as it is by the weak characterization.

Orthogonalization step

The orthogonalization step consists in searching the PC representation of Φ (k)⊥ r orthogonal to Φ q , when knowing the PC representation of Φ q and the one of Φ (k) r for q < r. For mass-normalized normal modes Φ q , the orthogonalization condition is defined by:

Φ T q M Φ (k)⊥ r = 0
and it is achieved for:

Φ (k)⊥ r = Φ (k) r -Φ q Φ T q M Φ (k) r ( 16
)
which is for the PC representation of the stochastic variables:

n λ φ=0 Ψ φ (ζ) φ (k)⊥ rφ = n λ φ=0 Ψ φ (ζ) φ (k) rφ - n λ φ=0 n λ ϕ=0 nm m=0 n λ λ=0 Ψ φ (ζ) Ψϕ (ζ) Ψm (ζ) Ψ λ (ζ) φ qφ φ T qϕ mmφ (k) rλ .
Then, by taking:

φ (k)⊥ ri = φ (k) ri -1 Ψ 2 i n λ ϕ=0 n λ φ=0 nm m=0 n λ λ=0 Ψ i Ψ φ Ψ ϕ Ψ m Ψ λ φ qφ φ T qϕ m m φ (k) rλ
for i = 0, . . . , n λ , it is ensured that:

n λ φ=0 Ψ i Ψ φ φ (k)⊥ rφ = n λ φ=0 Ψ i Ψ φ φ (k) rφ - n λ φ=0 n λ ϕ=0 nm m=0 n λ λ=0 Ψ i Ψ φ Ψ ϕ Ψ m Ψ λ φ qφ φ T qϕ m m φ (k) rλ .
However, the above expression for φ (k)⊥ rφ involves expectations of five PC products. Thus, an enhancedmore efficient-strategy is desirable from a numerical point of view to find φ (k)⊥ rφ . From:

n λ φ=0 Ψ φ (ζ i ) φ (k)⊥ rφ = Φ (k)⊥ r (ζ i ) ∀i
we know that:

φ (k)⊥ rφ = < Φ (k)⊥ r , Ψ φ (ζ) > < Ψ 2 φ (ζ) > or: φ (k)⊥ rφ = 1 < Ψ 2 φ (ζ) > ngp i=1 (Ψ φ (ζ i ) w i ) Ψ φ (ζ i ) φ (k)⊥ rφ (17)
when sampling points ζ i are the ones of the Gauss quadrature. To express

Ψ φ (ζ i ) φ (k)⊥ rφ
, it is known from the Eq. ( 16) that:

Φ (k)⊥ r (ζ) = n λ φ=0 Ψ φ (ζ) φ (k) rφ - n λ φ=0 n λ ϕ=0 nm m=0 n λ λ=0 Ψ φ (ζ) Ψϕ (ζ) Ψm (ζ) Ψ λ (ζ) φ qφ φ T qϕ mmφ (k) rλ (18)
and we state that the above equality holds if:

Ψ φ ζ j φ (k)⊥ rφ = Ψ φ ζ j φ (k) rφ - n λ ϕ=0 nm m=0 n λ λ=0 Ψ φ ζ j Ψϕ ζ j Ψm ζ j Ψ λ ζ j φ qφ φ T qϕ mmφ (k) rλ
for all ζ j . Hence, evaluating the above relation at sampling points ζ i enables the evaluation of φ (k)⊥ rφ from the Eq. ( 17).

For large structural systems, this strategy is suitable to construct orthogonalized statics or inertial modes that are necessary to enrich a truncated normal mode basis, as well as for modal synthesis methods [START_REF] Géradin | Mechanical vibrations: Theory and applications to structural dynamics[END_REF].

Modal projections

Modal projections of structural matrices are performed easily in a strong characterization. Considering a mass-normalized normal mode Φ q and the stiffness matrix K, the modal projection Λ q is defined as:

Λ q (ζ) = Φ T q (ζ) K (ζ) Φ q (ζ)
while its PC representation is:

Λ q (ζ) = n λ λ=0 Ψ λ (ζ) λ qλ for: λ qλ = < Λ q (ζ) , Ψ φ (ζ) > < Ψ 2 φ (ζ) > .
They are given numerically from a quadrature rule:

λ qλ = 1 < Ψ 2 λ (ζ) > ngp i=1 Λ q (ζ i ) Ψ φ (ζ i ) w i
where:

Λ q (ζ i ) = n λ φ=0 n k k=0 n λ λ=0 Ψ φ (ζ i ) Ψ m (ζ i ) Ψ λ (ζ i ) φ T qφ k k φ qλ .
An identical strategy enables the computation of modal damping.

Determination of the stochastic normal modes basis using the subspace inverse power method and polynomial chaos representation

Let us now take a look at the stochastic version of the subspace inverse power method presented in the previous subsection 4.1, to see if it is a suitable alternative to the deflated inverse power method in order to build the stochastic normal modes basis.

The proposed algorithm is the following. For an initial set of n r modes Φ (0) having the PC representation

Φ (0) , do: 1. Choose Φ RBM from null (k) to represent Φ RBM 2. Normalize Φ RBM such that Φ T RBM M Φ RBM = i 0 3. From φ q , define Φ (0)
q to represents Φ (0) q for normal modes corresponding to non null eigenvalues 4. Orthonormalize initial matrix Φ (0)

Evaluate Λ (0)

q from a modal projection of K over Φ (0) q 6. For k = 1, 2, . . . do:

(a) Let F (k) q = M Φ (k-1) q (b) Solve K Φ (k) q = F (k) q (c) Orthonormalize Φ (k) (d) Evaluate Λ (k) q from a modal projection of K over Φ (k) q 7. Let Φ q ← [Φ 0 , Φ q ] 8. End
In this algorithm, the sub-step (c) implies the orthogonalization of all normal modes, including the rigid body modes, when they exist. As for the stochastic deflated inverse power method, the starting guess Φ (0) is produced from the normal mode matrix φ of the deterministic nominal mechanical system obtained by the deterministic subspace inverse power method, and orthonormalization is achieved as in the above subsections 4.4 and 4.3:

1. Evaluate Φ (k)1 1
that normalizes modal mass to unity and let

Φ (k) 1 ← Φ (k)1 1 2. For r ← 2 to n r , (a) Find Φ (k)⊥ r that orthogonalizes Φ (k) r to Φ (k) p for p ∈ {1, . . . , r -1} and let Φ (k) r ← Φ (k)⊥ r (b) Find Φ (k)1 r
that normalizes the modal mass to unity and let

Φ (k) r ← Φ (k)1 r
3. End and the convergence condition for Λ q reads: max r∈{1, ..., nr}

Λ (k) r -Λ (k-1) r Λ (k) r < .
Comparisons of this algorithm with the one proposed for the stochastic version of the deflated inverse power method show that differences occur only in the organization of the numerical operations and resulting computation time but not for obtained results.

5 Numerical applications

First application

The first application concerns a system having two degree of freedom and involving three random variables.

Introducing ζ l for l = {1, 2, 3}, three uncorrelated standard normal variables, the random structural matrices are:

K = µ a 1 + ζ 2 1 -1 √ 2 σ a µ a    1 -1 -1 1   
and:

M =    1 0 0 1    + 1 √ 2 σ a µ a    ζ 2 2 -1 0 0 ζ 2 3 -1    .
The parameters are: µ a = 20 and σa µa = 5 %. Notice that positive definiteness is ensured for the stiffness and masses as long as σa µa < √ 2.

Hermite polynomials up to the second degree (10 terms) are used for the PC representation of the stochastic eigensolutions. Eigensolutions of the nominal system are:

φ = 1 √ 2    1 1 1 -1    and λ =    0 0 0 40   
showing that the first normal mode is a rigid body mode. Then, considering the stochastic eigenproblem, the first eigenvalue is deterministic, indeed it is null. The amplitude of the normalized rigid-body mode is however stochastic since it has a unit norm with respect to an inner product defined by the stochastic mass metric. The PC representation of this first normal mode is obtained by normalizing the starting guest which is built from the nominal system. The second normal mode is sought from the deflated inverse power method presented above, using the enhanced strategy for the orthogonal step. Since the iteration matrix K is not invertible, a Moore-Penrose pseudo inverse is involved for the generation of a new iterate. When starting Figure 1 shows the PDFs of the second eigenvalue and Figure 2 shows the PDFs of the eigenvectors components obtained. These results are compared to a MCS which has been carried out using 5 × 10 5 sample size. A satisfactory agreement is achieved.

Second application

The second application is inspired by the simplified three DOFs model of the bladed disk of the reference [START_REF] Ghosh | An invariant subspace-based approach to the random eigenvalue problem of systems with clustered spectrum[END_REF]. Introducing ζ l for l = {1, 2, 3, 4}, four uncorrelated standard normal variables, the random structural matrices are:

K = µ b       1 0 0 0 1 0 0 0 1       + σ b √ 2       ζ 2 1 -1 0 0 0 ζ 2 2 -1 0 0 0 ζ 2 3 -1       +µc 1 + ζ 2 4 -1 √ 2 σc µc       2 -1 0 -1 2 -1 0 -1 2      
and: The parameters used are the ones proposed in the reference [START_REF] Ghosh | An invariant subspace-based approach to the random eigenvalue problem of systems with clustered spectrum[END_REF]: µ b = 100, µ c = 20, while σ b µ b = σc µc = 5 % and Hermite polynomials are chosen for the PC representation of the stochastic eigensolutions and for the PC representation of the stochastic dynamic responses. The deflated inverse power method presented above is chosen to find the three stochastic normal modes while the normal modes of the nominal system are used as starting guesses. These normal modes are:

M =       1 0 0 0 1 0 0 0 1       + 1 √ 2 σ b µ b       ζ 2 1 -1 0 0 0 ζ 2 2 -1 0 0 0 ζ 2 3 -1       .
φ = 1 2       1 √ 2 1 √ 2 0 - √ 2 1 - √ 2 1       .
Reference results are produced by MCS carried out using 5×10 5 sample size and a Modal Assurance Criterion to identify and to order the stochastic normal modes. Empirical means and standard deviations of stochastic eigenvalues are given in Table 1 and their PDFs are plotted in Figure 3. Empirical CDFs of {Φ r } i , the components i of the normal modes Φ r , are shown in Figure 4. These MCS results are reference results to be compare with the one produced by the PC representation. For the current application, the 3rd component of the 3rd normal mode is the most critical regarding the convergence issue. For this normal mode component, if the 2nd degree (15 terms) or the 3rd degree (35 terms) is chosen for the polynomial order involved in PC representation, we get an identical relative error of 50% in the skewness and a relative error of 75% in the kurtosis of this normal mode component. But for a fourth degree (70 terms), we get a very satisfactory agreement between MCS and PC representation, as shown in Figure 4. In this case, the deflated inverse power algorithm presented above requires 10, 14, and 2 iterations to converge to a relative precision = 10 -6

for each stochastic eigenfrequency of the normal modes.

To compute frequency responses functions, the structural forces and damping matrices are

f (ω) =       1 0 0 0 1 0 0 0 1       and C =       1 0 0 0 1 0 0 0 1      
. Figure 5 and6 show i) the mean; ii) the 5 percentile; and iii) the 95 percentile of the amplitude of the nine frequency responses functions, as well as their second statistical moment. Notice that computing the mean amplitude requires all the terms of the PC representation, and not only the first one. Hence, agreement with MCS can only be achieved for an adequate truncation of the PC representation. For the stochastic eigenproblem, a satisfactory agreement is achieved. However, small discrepancies can be observed at ω 2 110 and ω 2 ∈ [165, 175]. This indicates that the order of truncation is not totally adequate in these ranges for the frequency responses, while it is adequate for the stochastic

Third application

The third application is the simplified three DOFs model of the bladed disk of the reference [START_REF] Ghosh | An invariant subspace-based approach to the random eigenvalue problem of systems with clustered spectrum[END_REF]. This application is interesting as two of the three deterministic normal modes have the same natural frequency.

Introducing ζ l for l = {1, 2, 3, 4}, four uncorrelated standard normal variables, the random structural matrices are:

K = µ b       1 0 0 0 1 0 0 0 1       + σ b √ 2       ζ 2 1 -1 0 0 0 ζ 2 2 -1 0 0 0 ζ 2 3 -1       +µc 1 + ζ 2 4 -1 √ 2 σc µc       2 -1 -1 -1 2 -1 -1 -1 2      
and M is the identity. The parameters are: µ b = 100, µ c = 20, while σ b µ b = σc µc = 5 % and Hermite polynomials up to the fourth degree are used for the PC representation of the stochastic eigensolutions and for the PC 5: Superimposition of frequency responses amplitudes of the three DoF system for the second application using the PC representation up to the 4th degree; Mean, 5% and 95 % lines of confidence interval are represented; From up to down is for the first to the last degree of freedom; From left-hand to right-hand is for an excitation applied from the first to the last degree of freedom; Thick grey lines are for MCS, dashed black lines are for PC representation Figure 6: Superimposition of second moment of frequency responses amplitudes of the three DoF system for the second application using the PC representation up to the 4th degree; From up to down is for the first to the last degree of freedom; From left-hand to right-hand is for an excitation applied from the first to the last degree of freedom; Thick grey lines are for MCS, dashed black lines are for PC representation

φ = 1 6       2 √ 3 √ 6 3 √ 2 2 √ 3 √ 6 -3 √ 2 2 √ 3 -2 √ 6 0       λ =       100 0 0 0 160 0 0 0 160      
.

showing that the second and the third natural frequencies are equals. When starting from these eigenvectors of the nominal system, the stochastic deflated inverse power method presented above find the three stochastic modes using 10, 82 and 11 iterations for a relative precision = 10 -6 for each eigenfrequency. Results of means and standard deviations of stochastic eigenvalues obtained are given in Table 2.

To compute frequency responses functions, the structural forces, and the damping matrices are

f (ω) =       1 0 0 0 1 0 0 0 1       and C =       1 0 0 0 1 0 0 0 1      
. Figures 7 and8 show i) the mean; ii) the 5 percentile; and iii) the 95 percentile of the amplitudes of the nine frequency response functions, as well as their second statistical moments. The obtained results are compared to a MCS which has been carried out using 5 × 10 5 sample size.

A satisfactory agreement is achieved, especially around the resonant frequencies, showing that eigenvalues multiplicity is handled effectively.

Fourth application

The fourth application is a realistic mechanical situation by considering a problem where more degrees of freedom than the number of sough normal modes while repeated eigenvalues exist for the nominal system.

Its purpose is to show the applicability of the strategies presented above to obtain modal based frequency responses in the low frequency range. This application addresses frequency response functions of a simplified tower guyed structure. More precisely the responses in the plane Oxy at the guys connection point P are of interest, see Fig. 9. Overall dimensions are such that eigenvalues are repeated and close to each other. For a guyed mast, wind and ice are the major loads which alter the nominal tension in guys. Thus, in this application, tensions in guys are Figure 7: Superimposition of frequency responses amplitudes of the three DoF system for the third application using the PC representation up to the 4th degree; Mean, 5% and 95 % lines of confidence interval are represented; From up to down is for the first to the last degree of freedom; From left-hand to right-hand is for an excitation applied from the first to the last degree of freedom; Thick grey lines are for MCS, dashed black lines are for PC representation Figure 8: Superimposition of second moment of frequency responses amplitudes of the three DoF system for the third application using the PC representation up to the 4th degree; From up to down is for the first to the last degree of freedom; From left-hand to right-hand is for an excitation applied from the first to the last degree of freedom; Thick grey lines are for MCS, dashed black lines are for PC representation supposed stochastic and we are interested with uncertainties quantification in frequency response functions of the P point.

The proposed simplified guyed mast structure consists of a vertical mast-column and 3 guys (Fig. 9).

The mast is pinned at its base and guyed in three directions at 120°angles in the plan. Guys are connected to the mast at the P point located at a height of h = 135 m. Distance between guy anchors and base tower is d = 145 m. Guys nominal tensions are t = 40, 000 N at guy anchors. Their Young's modulus is e = 165 GPa, mass density is ρ = 8, 100 kg/m 3 and area a = 3 × 10 -4 m 2 .

Guy is a tensioned cable and cable structures are nonlinear elastic systems with large displacements. If we consider the cable in the Oxz plane, the static equilibrium of the cable subject to tension and self weight leads to the profile of the cable which may be described by [START_REF] Irvine | Cable Structures[END_REF]:

z = t x q cosh q t x x + a 1 - t x q cosh (a 1 ) + h
where t x is the horizontal component of the cable tension, q = -9.81ρa is the distributed self weight of the cable and:

a 1 = arcsinh   qh 2t x sinh qd 2tx   - qd 2t x .
In addition, the tension at any point on the cable is given by:

t s = t x cosh q t x x + a 1 hence t = t x cosh q tx d + a 1 .
Once the static equilibrium of guys is achieved, the geometry of the system is well defined and additional displacements and strains due to dynamic loads are supposed to be small. Thus, a dynamic linear behavior is assumed. Guys are then modeled by a serie of 20 truss finite elements producing a mass matrix and a stiffness matrix composed by the sum of an elastic stiffness matrix and a geometric stiffness matrix [START_REF] Géradin | Mechanical vibrations: Theory and applications to structural dynamics[END_REF]. For the response in the Oxy plane, the mast contribution to the stiffness is neglected. However, partial mass of the mast act on node P. It is modeled as a single lumped mass, being 5, 750 kg. The resulting global assembly leads to a system having the form of Eq. ( 1) where the damping matrix is c = 0.05m and the number of active degrees of freedoms for the finite element model is n = 173. The deterministic modal analysis with nominal tensions exhibit repeated or close eigenvalues: λ 1 = λ 2 = 3.04 rd 2 /s 2 , λ 3 = 3.29 rd 2 /s 2 , λ 4 = λ 5 = 3.46 rd 2 /s 2 , see Fig. 10, associated to damped factors

1 2 √ λ1 ξ 1 = 1 2 √ λ2 ξ 2 = 1.43 %, 1 2 √ λ3 ξ 3 = 1.38 %, 1 2 √ λ4 ξ 4 = 1 2 √ λ5 ξ 5 = 1.34 %.
It is chosen to approximate the frequency responses functions by truncating arbitrary the modal basis to the five first normal modes.

For a tower guyed, wind and ice are the major loads. The wind can increase the tension in one guy while it can decrease the tension in another guy. In contrast, the ice will result in an increase for tensions in all guys. Then, the 3 tensions in guys are supposed to be random and correlation between them is a parameter of the study. For the 3 random variables, we choose arbitrary to consider a uniform law having a 4.5 % coefficient of variation (being the standard deviation over the mean) while the mean is fixed at the nominal value t. Then, three situations are considered for correlation coefficients of the random variables:

1. ρ ij = 0 for i = j, i ∈ {1, 2, 3} and j ∈ {1, 2, 3}.

2. ρ ij = -0.4 for i = j, i ∈ {1, 2, 3} and j ∈ {1, 2, 3}.

3. ρ ij = 1 for i ∈ {1, 2, 3} and j ∈ {1, 2, 3}.

Accordingly with the uniform variables, Legrendre polynomials are chosen for the PC representation. We chose to truncate the polynomial basis at an order 8, leading to 165 polynomials terms and coefficients and 28, 545 unknowns for the system of Eq. [START_REF] Ghanem | Efficient characterization of the random eigenvalue problem in a polynomial chaos decomposition[END_REF]. Figures 12,13 and 14 show the mean; and ii) the 5 and the 95 percentile of the amplitudes of the two horizontal frequency response functions, as well as their second statistical moments, for the correlation coefficients of cases 1 to 3 (respectively). The obtained results are compared to MCSs which have been carried out using 10 6 sample size. A satisfactory agreement is achieved, even around the resonant frequencies for all these 3 cases. Thus, meaningful results can be effectively achieved by the proposed strategies even for such complex situations which involve close and repeated eigenvalues.

Summary

This work concerns modal frequency responses of stochastic discretized linear dynamic problems. It is limited to symmetrical structural matrices and proportional damping. The PC representation is the chosen spectral approach for random quantities. Within this framework, from the normal basis construction up to the modal synthesis, we propose efficient formulations and strategies in order to maintain the computational effort similar to the one involved for stochastic linear static problems.

To treat the stochastic eigenproblem, instead of a newton-raphson strategy, adaptation of the deterministic deflated inverse power and subspace inverse power method is achieved in the stochastic framework. The Figure 9: Squetch of the guyed tower deterministic subspace inverse power method which deal effectively with close and repeated eigenvalues is chosen to initiate starting guess of both these iterative algorithms. Examination of involved steps of both these stochastic versions shows that they both have the advantage of having a main step which involves the resolution of a linear system similar to the one produces by the static case. Thus, as for linear stochastic static problems, required expectations are limited to three PC products for this computation step. Moreover, we have pay special attention to propose efficient numerical formulations or strategies for ensuring orthonormalization of the modal basis and for obtaining modal based frequency responses. They also use only expectations of PC products which are already used for stochastic static problems. All these propositions ensure the numerical efficiency of the stochastic modal based frequency response analysis.

Next, three simple applications and one realistic situation are studied. MCS is avoided to initiate iterations of the stochastic deflated inverse power algorithm. As a starting guesses, eigensolutions of the deterministic nominal problems are chosen. Through these applications, applicability of the proposed strategies to obtain modal based frequency responses is shown even for the realistic situation where there is more degrees of freedom than the number of sough normal modes while repeated eigenvalues exist for the nominal system. It is also shown that the proposed method can handle null eigenvalues. Comparisons of obtained results with the ones computed with MCS are satisfactory. In addition, it is found that frequency responses are more

  r

Figure 1 :

 1 Figure 1: Superimposition of the PDF of the second eigenvalue for the first application; thick grey lines are for MCS, dashed black lines are for PC representation

Figure 2 :

 2 Figure 2: Superimposition of PDFs of the two components (column-wise) of the two normal modes (row-wise) for the first application; thick grey lines are for MCS, dashed black lines are for PC representation

Table 1 :Figure 3 :

 13 Figure 3: Superimposition of PDFs of the three eigenvalues for the second application; thick grey lines are for MCS, dashed black lines are for PC representation

Figure 4 :

 4 Figure 4: Superimposition of CDFs of the three components (column-wise) of the three normal modes (rowwise) for the second application; thick grey lines are for MCS, dashed black lines are for PC representation

Figure

  Figure5: Superimposition of frequency responses amplitudes of the three DoF system for the second application using the PC representation up to the 4th degree; Mean, 5% and 95 % lines of confidence interval are represented; From up to down is for the first to the last degree of freedom; From left-hand to right-hand is for an excitation applied from the first to the last degree of freedom; Thick grey lines are for MCS, dashed black lines are for PC representation

Figure 12 :

 12 Figure 12: Case 1 results for the guyed tower of the fourth application using the PC representation up to the 8th degree (dash-dot black lines) vs. MCS (continuous thick grey lines); Up: mean of the amplitudes of the 2 frequency responses functions, as well as the 5% and 95 % lines of confidence interval; Down: second moment of the amplitude of the 2 frequency responses functions; Left subfigures concern the response in the x-direction due to an excitation in the x-direction; Right subfigures concern the response in the y-direction due to an excitation in the y-direction

Figure 13 :

 13 Figure 13: Case 2 results for the guyed tower of the fourth application using the PC representation up to the 8th degree (dash-dot black lines) vs. MCS (continuous thick grey lines); Up: mean of the amplitudes of the 2 frequency responses functions, as well as the 5% and 95 % lines of confidence interval; Down: second moment of the amplitude of the 2 frequency responses functions; Left subfigures concern the response in the x-direction due to an excitation in the x-direction; Right subfigures concern the response in the y-direction due to an excitation in the y-direction

Figure 14 :

 14 Figure 14: Case 3 results for the guyed tower of the fourth application using the PC representation up to the 8th degree (dash-dot black lines) vs. MCS (continuous thick grey lines); Up: mean of the amplitudes of the 2 frequency responses functions, as well as the 5% and 95 % lines of confidence interval; Down: second moment of the amplitude of the 2 frequency responses functions; Left subfigures concern the response in the x-direction due to an excitation in the x-direction; Right subfigures concern the response in the y-direction due to an excitation in the y-direction

Table 2 :

 2 Statistics of normal modes frequencies for the third application

	Normal mode number r	1	2	3
	Mean of frequency	99.75 159.4 160.8
	Standard deviation of frequency 2.51	4.55	5.54
	representation of the stochastic dynamic responses. Eigensolutions of the nominal system are:

However, care must be taken when natural frequencies are repeated. In this situation any combination of associated normal modes leads to another also valid normal mode; that is the associated eigenspace has dimension bigger than one. In such case, shapes of both modes of the previous example looks like a mix of the two original modes. Then, a specific procedure must be involved to avoid difficulties in following their evolution.

difficult to represent with PC than the stochastic normal modes.