domain: shs info

text:
Ariane Pinche

Thibault Clérice

Alix Chagué

Jean-Baptiste Camps

Malamatenia Vlachou-Efstathiou

Matthias Gille Levenson

Olivier Brisville-Fertin

Federico Boschetti

Franz Fischer

Michael Gervers

Agnès Boutreux

Avery Manton

Simon Gabay

Patricia O'connor

Wouter Haverals

Mike Kestemont

Caroline Vandyck

CATMuS -Medieval: Consistent Approaches to Transcribing ManuScripts A generalized set of guidelines and models for Latin scripts from Middle Ages (8 th -16 th )

Introduction

In 2016, the OCRopy tool [START_REF] Breuel | Ocropy: Python-based tools for document analysis and OCR[END_REF]) -designed as an open-source, easily trainable OCR engine-piqued the interest of digital humanists in deploying its application to manuscripts Camps 2017. While the majority of the community explored its usefulness for incunabula and fraktur texts [START_REF] Springmann | Ground Truth for training OCR engines on historical documents in German Fraktur and Early Modern Latin[END_REF], we delved into its potential beyond printed material. Between 2016 to 2019, the landscape of handwritten text recognition (HTR) underwent a profound transformation. The advent of general public GPUs for model training became feasible, and the availability of programming libraries like Tensorflow and PyTorch became standard. Notably, during this period, Kraken [START_REF] Kiessling | The Kraken OCR system[END_REF]) gained prominence, while Transkribus [START_REF] Kahle | Transkribus -A Service Platform for Transcription, Recognition and Retrieval of Historical Documents[END_REF]) emerged as a comprehensive solution, enabling manual annotation, model training and post-correction.

Along with the growing popularity of applications such as Transkribus and eScriptorium [START_REF] Stokes | The eScriptorium VRE for Manuscript Cultures -Classics@ Journal. en-US[END_REF], the use of HTR has seen significant growth in the digital humanities. Each project varies widely and is often inclined to tailor its own transcription methods based on specific objectives: from automatically transcribing cartularies [START_REF] Stutzmann | Handwritten Text Recognition, Keyword Indexing, and Plain Text Search in Medieval Manuscripts[END_REF] to extracting literary manuscript texts for further study [START_REF] Camps | Noisy medieval data, from digitized manuscript to stylometric analysis: Evaluating Paul Meyer's hagiographic hypothesis[END_REF][START_REF] Haverals | The Middle Dutch Manuscripts Surviving from the Carthusian Monastery of Herne (14th century): Constructing an Open Dataset of Digital Transcriptions[END_REF]. The need for a standardization of methods and datasets became evident in 2021 with the production of generic HTR models able to embrace diversified collections. The divergence in transcribing practices across projects prompted the creation of a set of unifying guidelines for transcribing Old French manuscripts from the 10 th to the 15 th century (Pinche 2022b). The collective dataset expanded through external additions, culminating in the development of the first generic model for Medieval manuscripts: the CREMMA medieval model (Clérice, Pinche, and Vlachou-Efstathiou 2022a). Subsequently, fuelled by multiple funding sources and collaborations, both the guidelines and models are extending their reach to encompass various languages, with a coming expansion for Old English and Middle Dutch.

The Consistent Approaches to Transcribing Manuscripts (CATMuS) guidelines, dataset, and model are the result of an international collaboration aimed at standardizing transcription practices for historic documents in Western Europe. CAT-MuS currently includes past and ongoing projects with collaborators from Europe and North America: France (CREMMA, HTRomance (Clérice, Chagué, et al. n.d.), Fabliaux (Pinche and Pierreville 2023), Liber [START_REF] Aruta | Liber[END_REF], Gallic(orpor)a, HTRogène1 ), Italy (HTRomance, HTRogene), Switzerland (Gallic(orpor)a, FoNDUE), Canada (DEEDS [START_REF] Gervers | DEEDS HTR Dataset[END_REF]). While CATMuS aims at merging transcription practices for HTR across documents from the Middle Ages to contemporary times, this paper specifically focus on the medieval and 15-16 th centuries early print (incunabula and print with gothic typefaces) from the datasets developed around the CATMuS Guidelines to help design general model(s) for HTR.

Guidelines

Transcribing involve a nuanced balance between fidelity to the source and interpretive choices. For texts in Anglo-Norman, for example, the historians' editions have neither accents nor apostrophes: we thus read Dangleterre and labbe, instead of modern d'Angleterre and l'abbé (Breuil 2019, Introduction, p. 5). On the contrary, French literary texts editions expect standardized punctuation, word segmentation, and accents. Moreover, while a human reader can handle the lack of homogeneity, it is a significant challenge for computers. The lack of standardization inflates character counts and introduces ambiguity. For example, the ‹p› with a stroke could be represented by the ‹ꝑ› (Unicode A751)or by a ‹ꝑ› (0554) from the Armenian alphabet. Such variations introduce noise during model training, leading to inconsistent predictions. Such issues have been found regarding a generic diachronic model for Spanish called Coloso español on Transkribus (Fradejas 2023) (cf. Figure 1).

In tandem with a seminar held at the École nationale des Chartesthe first guidelines for medieval French emerged, addressing documents spanning from the 10 th to 15 th centuries (Pinche 2022b). These guidelines aim to assist medievalists in translating manuscripts into a computer-readable format and in unifying their practices by: 1. restricting the characters used for specific purposes according to the MUFI 2 ; 2. applying a graphematic transcription3 ; 3. keeping abbreviations (and thus, reducing the part of the language specific traits the model has to learn).

Decisions on abbreviations, diacritics, and phonetic representations aim to minimize characters in the HTR model. In medieval sources, distinctions like ‹u›/‹v› and ‹i›/‹j› that are typically variations in form rather than distinct phonetic realizations are reduce to ‹u› and ‹i›.

Ongoing efforts include addressing new transcription requirements arising from diverse language domains, illustrated by challenges introduced by Latin in the CREMMA project and Spanish in the HTRomance project. The addition of Latin as added new difficulties as the representation of specific abbreviations for morphological suffix representations such as ‹ꝵ› for ‹rum› or ‹› for ‹bus›. The future addition of Middle English introduces the imperative need for new signs, like ‹ð› (eth) representing the voiced and voiceless dental fricative ([θ] and [ð]).

The establishment of transcription rules remains an ongoing and evolving effort, crucial for ensuring data quality and homogeneity. To ensure a better accessibility to the guidelines, the project aims at building a new set of guidelines, more interactive and more multilingual than the current PDF, with a delivery planned for the beginning of 2024.

Data

The CATMuS Medieval 1.0.0 dataset is built upon 17 different repositories across 10 different projects for which the guidelines have been applied, and for which Choco-Mufin character normalization and quality control have been applied Clérice and Pinche 2021. The corpus (c. 115.000 lines and 3.4M characters) is heavily influenced by its two first languages: Old French and Latin, and the existence of three personal projects4 which produces the three overrepresentation of French-Narratives (29.2% of the dataset), Treatises-Latin (24.2%) and Treatises-Castillan (23%, cf. Figure 2a). Speeches, Epistolary and Drama are three under-represented genres but it is not expected that their abbreviation rate or vocabulary differs much from the Narrative and Poetry genres.

Languages from the modern area of Italy and Spain (Castilian, Catalan, Navarrese and Venitian) are mostly underrepresented as only two projects are providing ground-truth. The final corpus contains 180 different documents-hands, of which 20 are prints and 34 are microfilms. Its most represented centuries are the 15 th century (74 document-hands), the 14 th (45), and the 13 th (34). The least represented centuries (8 th , 10 th , 11 th with a total of 6 manuscripts) can only be represented by Latin manuscripts and Old English manuscripts, which reduces the opportunities to accumulate data.

Evaluation

For CATMuS Medieval 1.0.0 [START_REF] Pinche | [END_REF], we provide a test based on out-of-domain manuscripts (see Table 1). Data were transcribed using the previous CREMMA Generic 1.0.1 model and post-corrected in the eScriptorium interface by the same annotators than the training data. Cambridge, CCC, ms. 041 represents one exception, as its language is completely out of domain (Old English and Middle Dutch is not in our training set), and its annotator has joined the initiative after the training of the model. Two competitor models are provided as baselines: the CREMMA Cortado mode (Pinche 2022a), a monolingual model trained on Old French (no character count available), and the CREMMA Generic model (Clérice, Pinche, and Vlachou-Efstathiou 2022b), trained on Latin and Old French (1.357 M characters).

Test results (cf. Figure 3) shows a large upgrade from the previous monolingual model (CREMMA Cortado) and a relative upgrade compared to the bilingual model CREMMA Generic. The later outperforms our larger models on most of the Old French test set (< 2%) except for Paris, BnF, fr. 1553. Looking at this specific manuscript, we hypothesize that its digitization (lower quality, a little blurry, greyscale microfilm) provided noise for our original model which our new model, trained with more noisy digitization, was able to overcome. On all other languages (except Old English and Middle Dutch), the model severely outperforms the original model.

For English specifically, we can see that moving from a monolingual to a multilingual model provides an important performance gain, while the language barrier remains stable for multilingual model, despite having more languages (+2) and more data (around 2 million more characters).

Overall, the new multilingual model provides a lower standard deviation across accuracies and a better average accuracy (Table 2), providing as such a better "generic" model according to this test set. The improvement from the Cortado model remains a huge leap, while the improvement from CREMMA Generic (and the detailed loss for some manuscripts) provides a picture of what can be enhanced in the near future.

Future

Consistent transcription standards, collaboration and data sharing across languages have provided HTR for medieval manuscripts with a boost in terms of accuracy. Such boost results in the provision of general models for HTR that play a crucial role in setting common benchmarks and improving the quality of predictions [START_REF] Camps | Make Love or War? Monitoring the Thematic Evolution of Medieval French Narratives[END_REF]. Thanks to the guidance provided by the guidelines and the support of our we have seen a significant increase in the number of available, interoperable datasets for training HTR models in recent years. This increase has enabled us to develop new generic models adapted to medieval manuscripts.

Although this article focuses mainly on Romance languages, we plan to extend the project to Old English and Middle Dutch, demonstrating our willingness to further broaden the linguistic eras treated. In addition, the project includes a section dedicated to the modern period, enabling the project not only to cover numerous linguistic eras, but also to extend chronologically from the Middle Ages to the contemporary era.

Latin

Figure 1 :

 1 Figure 1: Prediction using Transkribus Coloso Español model. The model uses two different ways to deal with abbreviations, line 24 ‹q› and line 27 ‹<ue>›. Capture kindly provided by J. M. Fradejas Rueda.

Figure 2 :

 2 Figure 2: Relative importance (in characters percentage) of categorical features in all CATMuS Medieval 1.0.0 training languages. Bars are accumulation of percentage per categorical feature.

Figure 3 :

 3 Figure 3: Test micro-accuracy per language and scripta. Bars represents the total number of characters per categorical feature.

Table 1 :

 1 Out-of-domain tests for various manuscripts in languages represented in the CATMuS initiative. As CATMuS 1.0.0 has no Old English manuscript nor Middle Dutch. It is expected that Cambridge, CCC, ms. 041 might have a drop in performances because of its specific characters: we provide a second scoring method ignoring thorn, eth and other Middle English specific characters. For Middle Dutch, confusion of u/n was a high factor in performance drop.

	Model

Table 2 :

 2 Scores across manuscripts depending on the model. Differences are computed against CATMuS. Micro average takes into account each manuscript character and error count, while Macro-average and STD take into account each manuscript accuracy.

This dataset is not yet public but collaborators were actively involved in the previous dataset.

https://mufi.info/q.php?p=mufi.

Transcription preserving letter sequences and reducing each form's meaning to an alphabetical system[START_REF] Stutzmann | Paléographie statistique pour décrire, identifier, dater[END_REF]

Namely the PhDs and MA thesis of Pinche (2021), Gille Levenson (2023b), and Vlachou-Efstathiou (2023).


Benoît Allignet

Camille Roukoz

Anne-Agathe Serre

Coralie Marjollet

Olivier Trédan

Pierre Meeus

Sébastien Mourey

Christine Rousset Jablonski

Isabelle Ray-Coquard

Disparités nationales de répartition des essais cliniques en onco-gynécologie

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

Introduction et objectifs :

Les femmes sont sous-représentées dans les essais cliniques en cancérologie [START_REF] Lee | Gender and sex disparity in cancer trials[END_REF][START_REF] Unger | Systematic review and metaanalysis of the magnitude of structural, clinical, and physician and patient barriers to cancer clinical trial participation[END_REF], alors même que moins de 5% des patients sont inclus [3]. Ces essais peuvent pourtant représenter des opportunités de traitements particulièrement prometteurs à l'échelle individuelle. Le principal déterminant de ce faible taux d'inclusion est l'indisponibilité des essais cliniques dans le centre de prise en charge [START_REF] Unger | Systematic review and metaanalysis of the magnitude of structural, clinical, and physician and patient barriers to cancer clinical trial participation[END_REF]. Nous avons ainsi analysé la répartition géographique des essais cliniques ouverts en France métropolitaine, en regard de la répartition médicale.

Méthode : Le nombre d'études-centre, défini par le nombre d'études multiplié par le nombre de centres investigateurs, a été recueilli à l'échelle départementale en France métropolitaine. Pour cela, une recherche a été effectuée le 30 août 2023 sur la base de données de l'application ScreenAct® (https://app.screenact.fr/), analysant la base ClinicalTrials.gov par traitement automatique du langage naturel. Le nombre d'études-centre a été rapporté par million d'habitants de sexe féminin (source: www.insee.fr/fr/statistiques/), pour les 3 importants groupes néoplasiques gynécologiques : cancers du col et/ou du vagin et/ou de la vulve, cancer de l'endomètre, et cancers tubo-ovariens. Etaient inclus tous les essais cliniques, observationnels et interventionnels, enregistrés sur ClinicalTrials.gov comme ouverts aux inclusions dans au moins un centre français. Les données démographiques médicales ont été obtenues sur l'atlas 2023 du Conseil National de l'Ordre des Médecins (source: www.conseil-national.medecin.fr/publications).

Résultat : En France métropolitaine, 37/96 départements (38.5%) ne disposent d'aucun essai clinique ouvert en onco-gynécologie. Les densités médianes d'études-centres ouvertes dans les cancers du col, de l'endomètre et de l'ovaire sont respectivement de 3.4, 1.5 et 2.3 par million d'habitants. La répartition en est en revanche très inégale (1er et 3ème quartiles, 0-14.4, 0-8.9, 0-11.4), la densité pouvant dépasser les 50 études-centres par million d'habitant (Tableau 1). Ces disparités sont d'autant plus marquées si on se focalise sur les essais interventionnels de phase 1 à 3 (Figure 1). Les répartitions géographiques des 3 groupes pathologiques sont superposables (Figure 1), hautement corrélés (Figure 2), mais ne reflètent que modérément la démographie médicale (Figure 2 et Figure 3Erreur ! Source du renvoi introuvable.). Discussion/Conclusion : La densité d'essai clinique en onco-gynécologie est extrêmement variable entre départements de France métropolitaine. Ceci peut limiter à la fois la représentativité des échantillons, conduire à des défauts d'inclusions et fermeture précoce d'essais, mais surtout limiter l'accessibilité à l'innovation médicale pour de nombreuses patientes habitant dans les déserts médico-scientifiques. La décentralisation de certains essais cliniques, l'organisation en réseaux d'excellences tel que FEM-NET, en réseaux de patients tels qu'IMAGYN, et l'utilisation de nouveaux outils numériques tels que ScreenAct® pourraient ainsi bénéficier à la fois aux patients et aux investigateurs.

Figure 1. Cartographies nationale métropolitaine de densité d'essais cliniques par million de femmes dans les différents groupes néoplasiques gynécologiques.

Abréviations : CHU, centre hospitalo-universitaire ; CLCC, centre de lutte contre le cancer

Figure 2 .

 2 Figure 2. Corrélogramme entre les études-centres des grandes pathologies onco-gynécologiques, les démographies médicales et la présence ou non de centres à vocation de recherche. Sont barrés les relations non significatives (p>0.05). Echelle d'intensité de couleur correspondant, à chaque intersection, au coefficient de corrélation (entre -1 et 1) entre la variable en abscisse et celle en ordonnée.Abréviations : Chir.Visc.Dig, chirurgien viscéral et digestif ; Gyn.Obs., gynécologue obstétricien ; Gyn.Med., gynécologue médical ; Onco.Med, oncologue médical ; Onco.RT., oncologue radiothérapeute ; CHU, centre hospitalo-universitaire ; CLCC, centre de lutte contre le cancer

Figure 3 .

 3 Figure 3. Cartographies nationale métropolitaine de densité de spécialistes (pour 100.000 habitants, tout sexe confondu) en (A) chirurgie viscérale et digestive, (B) gynécologie-obstétrique, (C) gynécologie médicale, (D) oncologie médicale et (E) oncologie radiothérapie, (F) anatomopathologie (d'après www.conseil-national.medecin.fr/publications)


Guilherme Viana
email: g.viana@avnir.fr

Guillaume Puel
email: guillaume.puel@centralesupelec.fr

Ludovic Chamoin
email: ludovic.chamoin@ens-paris-saclay.fr

Andrea Barbarulo
email: andrea.barbarulo@centralesupelec.fr

Implementation and analysis of viscoelastic damping in a 2D + 1D model of railway track vibrations

Keywords: Rail vibrations, SAFE method, Damping phenomena, Finite Element, Viscoelastic behavior

One of the most important technological challenges in the railway context nowadays is reducing the noise produced by the interaction between the train wheels and the track surface, that is the rolling noise. Introducing viscoelastic damping in the railway structure is one method to mitigate this rolling noise. In the paper, we investigate the impact of damping in the vibrational response of railway tracks using different viscoelastic models and parameters available in the literature. For this purpose, the main focus here is on developing an advanced numerical model in which viscoelastic models representing the rail foundation and/or pad are implemented in the framework of the Semi-Analytical Finite Element (SAFE) method. Linear or nonlinear viscoelastic models are considered, compared and effectively implemented for accelerance computation. The obtained numerical strategy leads to affordable and accurate predictions on dynamic features of the railway structure compared to classical methods.

Introduction

Reducing the noise level generated by railway traffic has been an important technological challenge of the last two decades, as well as an attractive research topic that gave rise to many scientific contributions (e.g. see [START_REF] Degrande | A numerical model for ground-borne vibrations from underground railway traffic based on a periodic finite element-boundary element formulation[END_REF][START_REF] Thompson | Railway Noise and Vibration -Mechanisms, modelling and means of control[END_REF] for an overview). It is still a crucial societal issue nowadays, due to the increase of rail transport in highly populated areas. A specificity of the railway noise is that it originates from various sources, encompassing aerodynamics phenomena as well as vibrations of the train structure, equipments and ground.

We here focus on the vibrations of railway tracks that entail rolling noise, which is the main noise source for train speeds between 80 km/h and 300 km/h. The dynamic behavior of the rails, in both vertical and lateral directions and over a large frequency range, is indeed a strong source of noise during the train passing through. Rail and wheel roughnesses generate vibrations of the rail and its environment (wheel, sleepers. . . ) and produce noise through acoustic radiation (Figure 1). In this context, modeling and simulating the vibrational response of the rail subjected to the dynamic train loading has received considerable interest and has progressively become a procedure of choice for optimal design purposes. Engineering models for predicting the rolling noise of conventional railways, including ballasted and slab tracks, are now well established [START_REF] Thompson | Railway Noise and Vibration -Mechanisms, modelling and means of control[END_REF]. Several approaches are available in the literature to analyze wave propagation along the rail, associated with analytical, semi-analytical, or numerical methods. They are based on fundamentals of structural dynamics under moving loads [START_REF] Fryba | Vibration of Solids and Structures Under Moving Loads[END_REF], and they usually consider the rail as an infinite domain along its longitudinal direction.

Figure 1. Rolling noise generation mechanism (from [START_REF] Thompson | Railway Noise and Vibration -Mechanisms, modelling and means of control[END_REF]).

On the one hand, many works have proposed analytical studies of the rail dynamical response by considering simple models based on beam theory, associated with various representations of the loading and foundation (e.g., using a continuous Winkler foundation model). They usually address the problem in the frequency domain with Fourier's transform [START_REF] Kenney | Steady state vibrations of beam on elastic subgrade for moving loads[END_REF][START_REF] Harrison | The dynamic response of railway track to high frequency vertical excitation[END_REF]. However, all these models have some drawbacks : (i) the coupling between torsion and horizontal bending vibrations is not taken into account; (ii) the deformation of the cross-section, which is not negligible for medium and high frequencies, is not considered; (iii) wave modes other than those predicted by the beam theory can not be observed (in particular for high frequency solicitations).

On the other hand, numerical full 3D analysis of the problem as performed in [START_REF] Heckl | Simulations and analyses of train-induced ground vibrations in finite element models[END_REF][START_REF] Wang | Free and forced vibration of repetitive structures[END_REF][START_REF] Wu | The effects of track non-linearity on wheel/rail impact[END_REF] may appear costly when considering the specific geometry of the studied waveguide-type rail structure, with uniform geometric properties of the cross-section along the longitudinal axis. It may also be computationally expensive when the frequency increases as the mesh size needs to be decreased, even when taking advantage of the periodic support, such as in the wave finite element method (WFEM). The WFEM method considers a single periodic element using a 3D FEM meshing [START_REF] Arlaud | Etienne Validation of a reduced model of railway track allowing long 3D dynamic calculation of train-track interaction Paris[END_REF]. Similarly to the WFEM method, the semi-analytical finite element (SAFE) takes into account the periodicity of the rail representing it with a single element. Nevertheless, in the SAFE method, the FEM mesh is applied to the constant cross section (2D) while analytical harmonic functions are considered in the longitudinal direction (1D).

In order to circumvent the previously mentioned drawbacks, we deal with the SAFE-based method in the present paper. It was originally developed in [START_REF] Gavric | Finite element computation of dispersion properties of thin-walled waveguides[END_REF][START_REF] Gry | Dynamic modelling of railway track based on wave propagation[END_REF] before being applied and enhanced in several contributions [START_REF] Sheng | Modelling ground vibration from railways using wave-number finite-and boundary-element methods[END_REF][START_REF] Loveday | Semi-analytical finite element analysis of elastic waveguides subjected to axial loads[END_REF][START_REF] Nilsson | A waveguide finite element and boundary element approach to calculating the sound radiated by railway and tram rails[END_REF][START_REF] Barbarulo | An Arnoldi reduction strategy applied to the semi-analytical finite element method to model railway track vibrations[END_REF]. It shares similarities with the so-called spectral finite element (SFE) approach developed in [START_REF] Finnveden | Evaluation of modal density and group velocity by a finite element method[END_REF]. In the SAFE method, the rail vibration is described by means of a waveguide model along the longitudinal direction z, while deformations of the 2D cross-section according to the (x,y) plane are captured by means of a finite element model. The resulting waveguide-FEM model thus yields a displacement field with a spatial evolution under the form φ(x, y).e jkz ; it is given analytically in the longitudinal direction (with harmonic description of waves travelling through the rail) and numerically in the transverse directions (for the description of section modes). Moreover, a periodic representation is given to the discrete rail support in the longitudinal direction, with simple mass-spring systems; the periodic structure of the track is further exploited using the Floquet principle [START_REF] Brillouin | Wave propagation in periodic structures[END_REF][START_REF] Mead | Wave propagation in continuous periodic structures: research contributions from Southampton, 1964-1995[END_REF], restricting the analysis to a portion (single periodic cell) of the rail with prescribed transfer matrix. Because of the constant cross-section of the rail, the applied FEM discretization of a 2D domain (in order to compute displacements in the cross-section) highly decreases the overall computational cost compared to full 3D numerical analyses. Nonetheless, the SAFE method leads to an accurate representation of vibration phenomena, and a quantitative analysis of the resulting noise, for a given configuration of the railway track infrastructure. Recently, the SAFE method was complemented with a specific numerical strategy based on an iterative second-order Arnoldi reduction (SOAR) technique [START_REF] Bai | SOAR: A second-order Arnoldi method for the solution of the quadratic eigenvalue problem[END_REF] in order to circumvent stability issues for low frequencies [START_REF] Barbarulo | An Arnoldi reduction strategy applied to the semi-analytical finite element method to model railway track vibrations[END_REF]. Here, we go one step further by extending the numerical approach to a richer class of material behaviors.

Railway track systems have been under continuous development and many different structural architectures have been introduced to fulfill specifications. In order to reduce railway track vibrations, elastic elements or additional masses are usually inserted in the vicinity of the rail (e.g., elastic rail pads) or under the track (e.g., non-ballasted floating slab track). They enable to reduce forces with excitation frequencies above the resonance frequency [START_REF] Wilson | Control of ground-borne noise and vibration[END_REF][START_REF] Nelson | Recent developments in ground-borne noise and vibration control[END_REF][START_REF] Ntotsios | A comparison of ground vibration due to ballasted and slab tracks[END_REF]. Nowadays, railway engineers aim at developing new means to decrease noise levels.

A promising technological advance is the use of elements made of dampening materials. Several recent works focus, indeed, to better understand, both theoretically and numerically, the role of viscoelastic damping in vibration control in railways applications [START_REF] He | Wanming An improved dynamic model of suspended monorail train-bridge system considering a tyre model with patch contact[END_REF][START_REF] Papangelo | Norbert Self-excited vibrations due to viscoelastic interactions[END_REF][START_REF] Nguyen | A new fuzzy-disturbance observer-enhanced sliding controller for vibration control of a train-car suspension with magneto-rheological dampers[END_REF][START_REF] Foyouzat | M In-plane vibration analysis of horizontally curved beams resting on visco-elastic foundation subjected to a moving mass[END_REF]. One possibility is to insert such materials in the rail support or foundation, replacing the traditional ballast. These can be modeled as linear or nonlinear viscoelastic dampers for the analysis. Railway track designers currently investigate how such elements could be advantageously used in order to minimize the acoustic pollution. Nevertheless, taking into account viscous effects with potential highly nonlinear behavior (which may significantly change the vibrational response of the railway track) in numerical modeling is often associated with the introduction of numerical complexities. Here again, most of the research works dealing with viscoelastic effects in the dynamical response of railway tracks use beam models. The dynamical analysis of beams resting on a viscoelastic foundation under moving loads has been extensively investigated. In [START_REF] Huang | Dynamic stiffness of infinite Timoshenko beam on viscoelastic foundation in moving coordinate[END_REF][START_REF] Huang | Response of an infinite Timoshenko beam on a viscoelastic foundation to a harmonic moving load[END_REF], the response of an infinite Timoshenko beam on a viscoelastic foundation to a harmonic moving load is studied by means of the dynamic stiffness matrix, as a function of the velocity and frequency of the load. In [START_REF] Pavlovic | Dynamic stability of Timoshenko beams on Pasternak viscoelastic foundation[END_REF], the dynamic stability of a Timoshenko beam supported by a Pasternak viscoelastic foundation subjected to compressive axial loading is investigated. Nevertheless, all these studies suffer from the same shortcomings listed above for models based on beam theory.

The main objective of this work is to effectively implement and analyze the capabilities of the SAFE numerical method in accounting for viscoelastic damping in the vibrational behavior of rail and associated accelerance, using viscoelastic parameters available in the literature. In order to get reliable information for comparison and optimization, we propose to extend the framework of the SAFE method in the presence of (possibly nonlinear) viscoelastic effects. We thus develop and implement a novel SAFE-based numerical approach similar to the method proposed in [START_REF] Barbarulo | An Arnoldi reduction strategy applied to the semi-analytical finite element method to model railway track vibrations[END_REF] and extend it to take into account a set of viscoelastic material models to represent the rail environment (support, foundation) in the semi-analytical formulation. This requires to adapt the solution scheme, in particular for the nonlinear case where a Newton-Raphson algorithm is additionally employed. The main change in the numerical scheme occurs when computing the impedance matrix of the rail support where a loop over the frequency domain is added in order to compute the complex stiffness of viscoelastic materials. Propagative wavenumbers and corresponding vibration modes of the rail are then computed as a function of the excitation frequency.

The outputs of the approach, in terms of vertical accelerances over the whole frequency range, permit to numerically and accurately quantify the action of dampening features of the rail environment. For illustration purpose, several viscoelastic models are considered even though the objective is not to compare them.

The paper is organized in five parts. After a literature review and introduction in Section 1, we introduce in Section 2 the railway structure and the SAFE-based model. In Section 3, we present the viscoelastic models used in this work. Next, we detail the implementation procedure of the viscoelastic material in the rail support, which is the originality of this paper. In Section 4, the performance of the approach is analyzed through several numerical experiments. Finally, conclusion and prospects of this work are addressed in Section 5.

Railway track model and Semi-Analytical Finite Element (SAFE) method

Railway track model

In this section we introduce the railway track model and its mechanical components. The railway track is made of two parts:

(i) a substructure, corresponding to the foundation, made of the ballast (for standard architectures) and possibly other layers (in the case of multi-layered foundation). The ballast may be replaced with concrete pavement in rigid tracks, which usually requires the additional use of viscoelastic elements in order to reduce noise emission;

(ii) a superstructure, called slab, located above the foundation. It is composed of the rail and the rail support, for instance the sleeper and rail pads.

The rail, in direct contact with the train wheels, has a specific constant transverse cross-sectional profile (see Figure 2). The geometry of this profile, as well as material properties of the rail, are standardized. The rail is attached to sleepers by means of rail pads and fastenings. These sleepers permit the transfer of forces between the rail and the foundation, as well as to maintain a constant distance between parallel rails. Rail pads and fastenings transfer and dampen the vertical load between the rail and sleepers. The railway track shown in Figure 2 has specific features : (i) the rail cross-section has constant physical and geometrical properties in the longitudinal direction, so that a waveguide model may be used to study wave propagation in the rail; (ii) the structure is periodic, so that a single periodic segment (cell) may be studied. In the following, we describe the SAFE numerical method that refers to such modeling ingredients.

Basic principles of the SAFE method

The SAFE method is based on a specific, cheap, but still accurate modeling of the wave propagation problem inside a waveguide-like rail with deformable cross-section as shown in [START_REF] Barbarulo | An Arnoldi reduction strategy applied to the semi-analytical finite element method to model railway track vibrations[END_REF][START_REF] Cettour-Janet | Modelling the vibrational response and acoustic radiation of the railway tracks PhD thesis from Université Paris-Saclay[END_REF]. The method developed here is similar to the one proposed in [START_REF] Barbarulo | An Arnoldi reduction strategy applied to the semi-analytical finite element method to model railway track vibrations[END_REF], including the use of the Arnoldi reduction strategy. The method focuses on a single periodic segment of the rail (Figure 2). The principle is to mix a FEM representation of the transverse deformation (i.e., over a 2D domain alone) and an analytical harmonic representation of the propagation along the longitudinal direction (taking advantage of the constant railway cross-section). In this section, the rail and its environment are supposed elastic; the viscoelastic behavior will be described in Section 3. The support with ballast is considered as a linear elastic component that translates and rotates, while rail pads are modeled by means of local equally spaced springs at the contact points along the rail foot and the sleeper (circular black points in Figure 3). To take into account the pad width, this latter is modelled with N n groups of springs. Those springs only account for translational displacements. The sleeper, linking the rail and the ballast, is supposed to be a rigid body. The steps of the SAFE method to compute the vibratory behavior of the rail are the following: (i) Computation of propagation modes in a free rail, that is without any loading or imposed boundary condition. After defining a FE mesh in the cross-section of the rail, free waves are determined from a quadratic eigenvalue problem (see below). Any displacement of the rail is then written as a linear combination of these wave modes; (ii) Introduction of boundary conditions by taking into account the structural rail environment (sleepers, rail pad and ballast), with linear elastic or viscoelastic coupling between the rail and the sleeper, and between the sleeper and the ballast. A transfer matrix is used to perform the coupling between adjacent periodic cells; (iii) Computation of the forced response with external loading applied to the rail.

Formulation of the SAFE method for the computation of free waves

We provide here basics on the SAFE formulation, but more detailed information about the derivation of equations of the SAFE method is found in [START_REF] Barbarulo | An Arnoldi reduction strategy applied to the semi-analytical finite element method to model railway track vibrations[END_REF][START_REF] Cettour-Janet | Modelling the vibrational response and acoustic radiation of the railway tracks PhD thesis from Université Paris-Saclay[END_REF]. The formulation of the SAFE method leans on the separation of variables between the transverse directions (x, y) and the longitudinal direction z. A displacement mode is written as: u(x, y, z, t) = φ(x, y).p(z, t)

in which p(z, t) = e j(kz-ωt) is the wave propagation function in the longitudinal direction, with k the wavenumber and ω the angular frequency. φ(x, y) is the 3D deformation mode shape of the cross-section, computed with FEM. It is written as φ(x, y) = N(x, y)Φ with N(x, y) the matrix constructed from 2D shape functions, and Φ the vector of nodal displacements.

Using the principles of the virtual works on a periodic element of length L, we obtain the dynamics equations of the rail :

B 0 p + [B 1 -B 2 ] ∂p ∂z -B 3 ∂ 2 p ∂z 2 + M ∂ 2 p ∂t 2 Φ = 0 B 2 p(0, t) + B 3 ∂p ∂z (0, t) Φ + G 0 - 1 2 F 0 = 0 -B 2 p(L, t) + B 3 ∂p ∂z (L, t) Φ + G L + 1 2 F L = 0 (2)
with M the mass matrix, G 0/L and F 0/L vectors associated respectively to the force reaction by the sleeper and the force from adjacent elements at the extremities. The last two equations from (2) correspond to force projections at both ends of the rail section. With p(z, t) = e j(kz-ωt) (straight track) and κ = jk, the first relation reads:

B 0 -ω 2 M + κ[B 1 -B 2 ] -κ 2 B 3 Φ = 0 or K 0 + κK 1 + κ 2 K 2 Φ = 0 (3)
with K 0 = B 0 -ω 2 M and K 2 = -B 3 some symmetric matrices (in the undamped case), and K 1 = B 1 -B 2 a skew-symmetric matrix. Equation (3) drives wave propagation in a free rail. It is a quadratic eigenvalue problem whose solutions represent free wave modes in an infinite rail free of external actions and support [START_REF] Gry | Dynamic modelling of railway track based on wave propagation[END_REF].

The couples (k, Φ) solutions of the system (3) may be computed for each value of the angular frequency ω, making calculations expensive. An alternative model reduction strategy is used in [START_REF] Gry | Dynamic modelling of railway track based on wave propagation[END_REF][START_REF] Barbarulo | An Arnoldi reduction strategy applied to the semi-analytical finite element method to model railway track vibrations[END_REF], requiring computations for a single arbitrary frequency value (e.g. 100 Hz). The N couples (k, Φ) associated with least decaying waves are selected to form a basis in order to compute solutions for other frequencies, with projection of matrices K j over the defined basis. In [START_REF] Gry | Dynamic modelling of railway track based on wave propagation[END_REF], it is shown that using modes computed at a 100 Hz frequency leads to an accurate representation of the track vibration phenomena over a frequency interval [50,5000] Hz. In [START_REF] Cettour-Janet | Modelling the vibrational response and acoustic radiation of the railway tracks PhD thesis from Université Paris-Saclay[END_REF], the basis is computed using the second-order Arnoldi reduction. In this case, the reduced basis is independent of the frequency value. The same strategy is used in the present work.

The general form of the displacement of the rail cross-section in the longitudinal position z is thus written as:

U(z) = N j=1 (α j e kj z ϕ j + β j e -kj z ψ j ) = F ϕ E(z)W α + F ψ E -1 (z)W β (4)
with F ϕ and F ψ some matrices whose columns are ϕ j and ψ j , representing the separation of waves propagating in positive and negative directions. The matrices W α , W β and E(z) are diagonal matrices built with α j , β j and e kj z , respectively. Coefficients α j and β j are obtained from a global computation taking into account boundary conditions and forces applied on the rail.

Computation of free waves in an infinite rail with periodic supports

The periodic feature of the support (represented by sleepers and rail pads) is taken into account in the SAFE method. As shown in Figure 3, the sleeper is represented by a rigid mass. The connection between the sleeper and the ballast considers the 6 degrees of freedom of the sleeper (3 translations + 3 rotations) whereas the connection between the sleeper and the rail only permits translations in the 3 directions. The force of the sleeper on the rail is obtained from the following steps:

(i) the dynamic equations (translation and rotation) for the sleeper are written:

-K b U + Nn i=1 K p i (U r i -U s i ) = -m s ω 2 U ; -C b θ + Nn i=1 P i K p i (U r i -U s i ) = -ω 2 J s θ (5)
with K b the translation stiffness matrix of the ballast, K p i the translation stiffness matrix of a pad, m s the sleeper mass, C b the rotation stiffness matrix of the ballast, and J s the sleeper matrix of inertia. Vectors U and θ represent displacements and rotations of the sleeper, respectively. Vectors U r i and U s i represent displacements of contact points of the rail and sleeper, respectively. The matrix P i is the rotational operator at contact points of coordinates (x i , y i , z i ); the link between the displacement of contact points of the sleeper and the generalized displacement/rotation of the sleeper is U s i = U -P i θ. (ii) from the previous relations, and assuming that K p i = K p /N n (i = 1, . . . , n), the displacement and rotation of the sleeper are calculated:

U = Nn i=1 Z U i K p U r i ; θ = Nn i=1 Z θ i K p U r i (6)
with

Z U i = A -1 (I -BZ θ i ) ; Z θ i = (D -CA -1 B) -1 (P i -CA -1 ) A = K b + N n K p -m s ω 2 I ; B = C T = -K p ( Nn j=1 P j ) ; D = C b + Nn j=1 P j K p P j -ω 2 J s (7)
(iii) from the displacement and rotation of the sleeper, the force generated by pad spring i on the sleeper is recovered as:

f i = K p (U r i -U s i ) = Nn j=1 K p (δ ij I -Z U j K p + P i Z θ j K p )U r j (8)
(iv) the global force F(z) of the sleeper on the rail is eventually computed as F(z) = A N T f (z)dA. It is related to the rail displacement through the impedance matrix of the rail support A p , such that

F(z) = A p U.
Then, in order to compute the free waves for the rail with periodic support, two ingredients are employed:

• The continuity of the displacement and equilibrium of forces at the interfaces between periodic cells Ω n and Ω n+1 are written:

U (n+1) (0) = U (n) (L) ; F (n+1) 0 = -F (n) L (9)
Using ( 2), (4), and the matrix operator A p , it yields:

F ϕ W (n+1) α + F ψ W (n+1) β -F ϕ E(L)W (n) α + F ψ E -1 (L)W (n) β = 0 F ϕ D k W (n+1) α -F ψ D k W (n+1) β -F ϕ D k E(L)W (n) α -F ψ D k E -1 (L)W (n) β -B -1 3 A p F ϕ E(L)W (n) α + F ψ E -1 (L)W (n) β = 0 (10)
where D k is a diagonal matrix containing wavenumbers k j . These relations lead to the following transfer equation:

γ (n+1) = (I + H -1 G)Lγ (n) = Tγ (n) (11)
with

γ (n) =   W (n) α W (n) β   ; H =   F ϕ D k -F ψ D k F ϕ F ψ   ; G =   B -1 3 A p F ϕ B -1 3 A p F ψ 0 0   ; L =   E(L) 0 0 E -1 (L)   (12)
T = (I + H -1 G)L is the transfer matrix between two elements.

• the Floquet principle γ (n+1) = λγ (n) , with scalar λ depending on frequency, is again used to relate solutions between adjacent periodic cells. It yields the eigenvalue problem:

(I + H -1 G)Lγ = λγ (13)
which may again be seen as a quadratic eigenvalue problem, and whose solutions permit to determine wave propagation features along the rail. The solution for a periodic supported rail with applied force is detailed in [START_REF] Barbarulo | An Arnoldi reduction strategy applied to the semi-analytical finite element method to model railway track vibrations[END_REF][START_REF] Cettour-Janet | Modelling the vibrational response and acoustic radiation of the railway tracks PhD thesis from Université Paris-Saclay[END_REF].

Introduction of viscoelastic behaviors in the SAFE method

In this section, we present the viscoelastic models used in this work and the implementation in the SAFE method.

Models of viscoelastic dampers

Viscoelastic dampers are widely used in engineering applications in order to control vibration of structures and associated noise [START_REF] Nashif | Vibration Damping[END_REF]. In this section, we shortly present the viscoelastic models considered in this work and their parameters. In order to implement a viscoelastic damper in the SAFE method, one needs to derive its frequency-dependent complex stiffness. We note E = E r (1 + jη) the complex stiffness, whose real part E r is the storage modulus, η is the loss factor and the product E r η is the loss modulus. Here, we consider that the translation and rotation stiffnesses have the same numerical value (but not the same dimension). We also choose an isotropic behavior which means that the stiffness in the viscoelastic models is independent of the direction. Table 1 shows the complex stiffness of the following models : Zener, generalized Maxwell model (GMM), Golla-Hughes-McTavish (GHM) and fractional derivative (FD) [START_REF] Baz | Active and passive vibration damping[END_REF]. The parameters needed to compute the complex stiffness in Table 1 are: k-stiffness, ω n -natural frequency, ρ-relaxation time, ζ-damping ratio, α-order of derivation (FD model), c-dissipation coefficient, andα n -gain (GHM model).

Another model implemented in this work is the cubic displacement stiffness (Duffing oscillator, we refer to this model as nonlinear (NL)), whose associated behavior is written as

F (t) = k L u(t)+c L u(t)+k N L u 3 (t)
where the parameter k N L is the stiffness related to the nonlinear term. This type of model is used in several

Model

Complex stiffness

Zener E ze = k 1 + k 2 jωρ 2 1+jωρ 2 GMM E GM M = k 0 + n m=1 kn jωρm 1+jωρm GHM E GHM = k 0 1 + N n=1 αn ω 2 +2ζnωnω ω 2 +2ζnωnω+ω 2 n FD E F D = k 0 +k∞(ρjw) α 1+(ρjw) α
works in which the ballast is modeled as a viscoelastic foundation. Dealing with a cubic displacement stiffness in the model requires the use of numerical methods for nonlinear solution such as Picard or Newton-Raphson methods (see section 3.3).

Viscoelastic models parameters

Some references determined the dynamic railpad and ballast properties for materials applied in the railway context. For example, the parameters of a Kelvin-Voigt model are estimated between 20 Hz and 2500 Hz in [START_REF] Maes | Measurements of the dynamic railpad properties[END_REF] . The stiffness and damping of elastomeric pads and ballast, in the frequency 20 Hz and 2500 Hz, are derived in [START_REF] Thompson | The dynamic behavior of rail fasteners at high frequencies[END_REF]. Due to the lack of data concerning viscoelastic materials applied in the railway structure in the frequency range of our interest (50-5000 Hz), in this work we adopted parameters found in the literature concerning viscoelastic materials applied in the engineering context for damping of structures, notably from [START_REF] Baz | Active and passive vibration damping[END_REF]. Those parameters are derived from experimental data within the frequency range considered here and they are in the same order of magnitude as the values reported in [START_REF] Maes | Measurements of the dynamic railpad properties[END_REF][START_REF] Thompson | The dynamic behavior of rail fasteners at high frequencies[END_REF].

The parameters for the FD model, α = 0.7, ρ = 0.004 s, k 0 = 3.62×10 5 N/m , k ∞ = 1.24×10 8 N/m are adopted from [START_REF] Baz | Active and passive vibration damping[END_REF]. The parameters for the Duffing oscillator are adopted from experimental data in [START_REF] Kargarnovin | Response of beams on nonlinear viscoelastic foundations to harmonic moving loads[END_REF]: In the case of GHM and GMM models, the parameters (properties and number of oscillators/elements) are computed from different sets of experimental data with different viscoelastic materials, given in [START_REF] Baz | Active and passive vibration damping[END_REF] on complex stiffness (storage modulus and loss factor from tests performed in the frequency range of interest in this work (50-5000 Hz)). We use experimental data for viscoelastic materials to estimate parameters using a least-square minimization. For the GMM model, the number of Maxwell elements in parallel for the best minimization of the least-square objective function is 6. The obtained parameters are shown in Table 2. The relaxation time values were chosen experimentally from [START_REF] Baz | Active and passive vibration damping[END_REF] so that to cover the whole measurement time.

k L =
The calculated parameters for the GHM are shown in Table 3. Three oscillators are necessary to obtain the best minimization of the least-square objective function [START_REF] Renaud | Rémi A new identification method of viscoelastic behavior: Application to the generalized Maxwell model Mechanical Systems and Signal Processing[END_REF]. In this minimization we adopted ζ n = 1 [START_REF] Baz | Active and passive vibration damping[END_REF]. Concerning the Zener model, we used the parameters values from reference [START_REF] Li | Andrea Testing and modelling of elastomeric element for an embedded rail system Materials[END_REF], that are k 1 = 62×10 6 N/m, k 2 = 9.6 × 10 6 N/m and c 2 = 7 × 10 5 Ns/m for its dissipation coefficient, with ρ 2 = c 2 /k 2 yielding 0.073.

The storage modulus and loss factor obtained from equations in Table 1 are shown in Figure 4. We observe that the storage modulus is almost constant except for the GHM model. The GMM model has the highest storage modulus value. As regards the loss factor, the Zener model has a constant and almost null value, whereas we observe a decrease of loss factor with respect to frequency for the FD and GMM models. The loss factor of the GHM models increases considerably from 50 Hz to 3000 Hz and then decreases slightly at 5000 Hz.

Implementation of viscoelastic models in rail support

The novelty in the present work lies in applying viscoelastic damping in the mechanical components of the rail such as the ballast and pad in the SAFE method. The use of a viscoelastic model implies the modification of the impedance matrix of the rail A p introduced in Section 2.4. In the case one wants to introduce one of the linear viscoelastic models or the nonlinear FD model in the ballast, the respective viscoelastic complex stiffnesses are introduced directly in (5), with K b = E I et C b = E I, with E being the rotational complex stiffness. We consider that the translational complex stiffness E and the rotational complex stiffness E have the same numerical value. In this work, for the pad we use the linear elastic behavior with reference parameters from Table 4 or the Kelvin-Voigt model as proposed for instance in [START_REF] Ntotsios | A comparison of ground vibration due to ballasted and slab tracks[END_REF] with K p = E KV I. The identity matrix renders the hypothesis that the translation stiffness is similar in all translation directions and the rotation stiffness is similar for all rotation directions in the ballast. In the case of the Duffing oscillator considered in the ballast and the Kelvin Voigt model in the pad, the dynamic equations for the sleeper (see equation ( 5)) are modified and the system of nonlinear equations of the sleeper's displacement and rotation as function of the rail displacement reads:

m s üx + (k L + N n k p x ) u x + c L ux + Nn i=1 k p x (-y i θ z ) + k N L u 3 x = Nn i=1 k p x u r xi m s üy + k L + N n k p y u y + c L uy + Nn i=1 k p y (x i θ z ) + k N L u 3 y = Nn i=1 k p y u r yi m s üz + (k L + N n k p z ) u z + c L uz + Nn i=1 k p z (y i θ x -x i θ y ) + k N L u 3 z = Nn i=1 k p z u r zi J x θx + e L + Nn i=1 y 2 i k p z θ x + d L θx + Nn i=1 k p z y i (u z -x i θ y ) + c N L θ 3 x = Nn i=1 k p z y i u r zi J y θy + e L + Nn i=1 x 2 i k p z θ y + d L θy - Nn i=1 k p z x i (u z + y i θ x ) + c N L θ 3 y = Nn i=1 -k p z x i u r zi J z θz + e L + Nn i=1 k p y x 2 i + k p x y 2 i θ z + d L θz + Nn i=1 k p y x i u y -k p x y i u x + c N L θ 3 z = Nn i=1 k p y x i u r yi -k p x y i u r xi . (14)
with k p z = k p y = k p x the direction components of the complex stiffness in the pad, u i and θ i being respectively the displacement and rotation in each direction (i = x, y, z). The parameters k L et e L are respectively the translation and rotational stiffnesses and they have the same numerical value, e L = 35.03 × 10 6 N.m/rad (see section 3.2 for k L ), even though not the same dimension. Moreover, d L is the rotational damping coefficient and c L the translational damping coefficient. In addition, we consider similar numerical values for those coefficients, with d L = 1732.50 × 10 3 N.m.s/rad (see section 3.2 for c L ). We also consider that the translational and rotational stiffnesses in the nonlinear terms in the Duffing oscillator have the same numerical value, that is c N L = 4.10 14 N.m/rad 3 (see section 3.2 for k N L ). Among the various methods to solve nonlinear equations such as [START_REF] Nilsson | A waveguide finite element and boundary element approach to calculating the sound radiated by railway and tram rails[END_REF], the force harmonic balance method [START_REF] Kovacic | The Duffing Equation: Nonlinear Oscillators and their behavior[END_REF][START_REF] Worden | Nonlinearity in Structural Dynamics[END_REF] is the most popular when the frequency domain is considered. The sleeper displacement solution of the nonlinear equation is written in the following form:

u(t) = N k=1 (a k cos(kωt) + b k sin(kωt)) . ( 15
)
The sleeper displacement in the complex domain is written as u(t) = U e jωt with U = a 2 + -bj 2 . In order to solve [START_REF] Nilsson | A waveguide finite element and boundary element approach to calculating the sound radiated by railway and tram rails[END_REF], we can apply the harmonic balance method with a single term for each component of displacement and rotation: u x = a 1 cos(ωt) + a 2 sin(ωt); u y = a 3 cos(ωt) + a 4 sin(ωt) u z = a 5 cos(ωt) + a 6 sin(ωt); θ x = a 7 cos(ωt) + a 8 sin(ωt) θ y = a 9 cos(ωt) + a 10 sin(ωt); θ z = a 11 cos(ωt) + a 12 sin(ωt). [START_REF] Finnveden | Evaluation of modal density and group velocity by a finite element method[END_REF] The nonlinear system J (a) = 0 obtained with ( 14) and ( 16) has 12 unknowns assembled in a vector a = [a 1 ...a 12 ]. For each degree of freedom (translation and rotation) there are two constants to calculate and the Newton-Raphson method is applied. For a given frequency ω i , the iterative method shown in Algorithm 1 depends on the initial solution a 0 that is taken as the solution of the previous converged vector (a i-1 ). The Jacobian matrix is noted J. Algorithm 1 Computing A p (ω i ) using a viscoelastic model for the (i)-th frequency Input: k = 0 and a 0 (a 0 = 0 for i = 1)

Compute K p = E KV (ω i )I or K p = E L for the pad if GMM, GHM, Zener or FD model then K b = E (ω i )I and C b = E (ω i )I, then solve (5) to compute U(ω i ) else (Duffing oscillator, see [START_REF] Nilsson | A waveguide finite element and boundary element approach to calculating the sound radiated by railway and tram rails[END_REF])

while e > 10 -6 do

k = k + 1 a (k+1) = a (k) -J(a (k) , ω i ) -1 J (a (k) , ω i ) e = a (k+1) -a (k) a (k) end while Compute U(ω i ) a 0 = a (k+1) initial guess for computing A p (ω i+1 ) for i > 1 end if solve F(z) = A p (ω i )U(ω i ) to compute A p (ω i )
In the numerical implementation, with respect to the linear elastic behavior, the introduction of a viscoelastic model requires adding a loop in frequency over Algorithm 1 in order to compute the frequencydependent complex stiffness.

Numerical results

In the first part of this section, the parameters used in the SAFE numerical model are introduced. Secondly, an analysis of the results obtained with the SAFE method in the linear elastic case (reference case) for the rail dynamic behavior is conducted for its mid-span cross-section and above the sleeper crosssection. We analyze the frequency behavior using as measure the accelerance, defined as the acceleration per unit of input force. In the third part of this section, the frequency behavior of a sleeper-viscoelastic ballast system is computed in order to investigate how each viscoelastic model behaves in a 1D vertical vibration problem. The properties of the viscoelastic models are those found from the minimization process in Section 3.2 or directly from the literature. Then, we apply some viscoelastic models only in the ballast using the SAFE method and we compare with a Timoshenko beam model. Next, we analyze the vertical accelerance response by introducing a viscoelastic behavior only in the ballast and then in the ballast and in the pad (using a Kelvin-Voigt model for the pad). Finally, a general analysis of accelerance variation is performed for the entire longitudinal length of a SAFE element by introducing the viscoelastic models in the ballast and in the pad. The number of modes for the reduced basis and the number of elements in the cross-section mesh (for type UI60 rail) are 50 and 137, respectively. Those values guarantee modes and mesh convergence for elastic and viscoelastic models.

Reference linear elastic properties and parameters in the SAFE model

The properties of the SAFE components and linear elastic ballast obtained from [START_REF] Cettour-Janet | Modelling the vibrational response and acoustic radiation of the railway tracks PhD thesis from Université Paris-Saclay[END_REF] and used as the linear elastic reference solution are described in this section. Table 4 shows the linear elastic properties of the ballast and pad. Table 5 shows the rail and sleeper properties.

Table 4

Linear elastic ballast and pad properties [START_REF] Cettour-Janet | Modelling the vibrational response and acoustic radiation of the railway tracks PhD thesis from Université Paris-Saclay[END_REF] Direction K b (N/m) C b (N.m/rad) loss factor : K b and C b K p (N/m) loss factor : K p x 10 Rail and sleeper contact points are the nodes at the rail foot for a given mesh. In this work, the crosssection mesh has 137 elements and 15 nodes at rail foot, as shown in Figure 5. The point location at which an effort is applied and the response is measured is the top of the rail. This point represents the contact between the rail and the wheel of a train (representing the rolling noise source).

Rail frequency response analysis using linear elastic ballast and pad given by SAFE method

In this section we analyze the frequency behavior of the the top of the rail cross-section (see Figure 5). Figure 6 shows the vertical accelerance measured at the top of the rail at mid-span and above the sleeper of a SAFE element using the elastic properties and parameters from Tables 4 and5. Four frequency regions limited by remarkable frequency peaks can be distinguished :

(i) In the first region, below the first frequency peak, the rail response is controlled by the ballast stiffness. The first frequency peak represents the rail and the ballast oscillating in phase without rail cross-section deformation. In that region, the rail response is similar at mid-span and above the sleeper. (ii) In the second region, between the first and second frequency peaks, the rail response is governed by the sleeper. There is an anti-resonance frequency where the accelerance decreases substantially. This corresponds to the system ballast-sleeper-pad acting as an absorber. The second frequency peak corresponds to the sleeper and rail vibrating in opposition of phase. (iii) In the third frequency region, the rail response is progressively uncoupled from the sleeper and ballast. The third frequency peak represents the 'pinned-pinned' frequency, where the semi-wave length of bending waves is the same as the distance between the sleepers. Above the sleeper, the pinned-pinned point is an anti-resonance mode. (iv) Finally, beyond the third frequency peak, the response of the rail is free from the sleeper and ballast influence. The cross-section deformation and pad stiffness govern the fourth region.

Study of a sleeper-viscoelastic ballast system

In order to study the vertical frequency response of the viscoelastic models, a simple 1D sleeperviscoelastic ballast system is proposed. The sleeper mass is given in table 5. Figure 7 shows the vertical accelerance response of the sleeper in a frequency range 50 Hz -5000 Hz. The resonance frequencies of the Zener, FD and GHM models are, respectively, about 150 Hz, 200 Hz and 700 Hz. We observe that the sleeper response using the Duffing oscillator (NL model) ballast and the GHM ballast are close to linear. The different responses observed are due to distinct complex moduli depicted in Figure 4. At 50 Hz, the accelerance depends on the complex moduli value. For example, the accelerance is higher for the GHM model as its stiffness is smaller compared to the other viscoelastic models. As the frequency is increased, we can observe the effect of dissipation (loss factor effect). For example, for the GHM model, the loss factor increases the damping and there is no peak in the accelerance. Nevertheless, for the GMM and FD models, the reduction in the loss factor with respect to frequency fosters the accelerance peaks. We consider the infinite Timoshenko beam model with periodic support proposed in [START_REF] Thompson | Railway Noise and Vibration -Mechanisms, modelling and means of control[END_REF] in order to compare with the SAFE model. The distance between two supports is 0.6 m, the same as what is adopted for the SAFE simulations. Figure 8 shows the comparison for the beam model with the SAFE model at mid-span of a SAFE element considering four different viscoelastic models. For each one, we observe the four different frequency regions, as explained in section 4.2. Below the pinned-pinned frequency, the vertical accelerance is similar between the beam and SAFE models. Above the pinned-pinned frequency, the accelerances peaks are different. This is explained by the deformation of the rail cross-section taken into account in the SAFE model [START_REF] Barbarulo | An Arnoldi reduction strategy applied to the semi-analytical finite element method to model railway track vibrations[END_REF], showing its relevance for the analysis of rail vibration in the frequency range considered here. The cross-section deformation appears as a deformation of the rail foot in a vertical movement (called "foot flapping", foot bending as a cantilever) or as a deformation of the rail web in the case of a lateral movement ("web bending") [START_REF] Barbarulo | An Arnoldi reduction strategy applied to the semi-analytical finite element method to model railway track vibrations[END_REF].

Viscoelastic models applied only in the ballast

The rail vertical accelerances using viscoelastic models only in the ballast at mid-span of the SAFE element are shown in Figure 9 (a). A significative vertical response of the rail due to viscoelastic ballast effect is observed particularly in the low-frequency range. The first frequency peak corresponding to rail and sleeper vibration in phase is dictated by the system ballast-sleeper shown in Figure 6 (a). The first anti-resonance frequency is observed for all models except for the nonlinear one (NL) 1 . For the GMM model, this point is quite shifted from the other models and occurs around 700 Hz. The second frequency peak corresponding to sleeper and rail vibrating in opposition of phase is similar for Zener, GHM, FD and linear elastic models. From a frequency of about 1000 Hz to 5000 Hz, the vibration behavior is similar for elastic and viscoelastic models and the pinned-pinned frequency is similar. As explained in Figure 6(a) for the sleeper-ballast model, the influence of the frequency dependent complex moduli is also observed in Figure 9(a), notably at low frequency, corresponding to the first region depicted in Figure 6(b). For example, the accelerance at the first frequency peak for the GHM model is less pronounced compared to the other models due to its important dissipation (increasing loss factor with respect to frequency, see Figure 4(b)). On the other hand, for the Zener model, the high accelerance peak at the first frequency peak is due to its low dissipation, as depicted in Figure 4(b). A similar pattern is observed for the vertical accelerance above the sleeper (Figure 9(b)).

The vertical accelerance in a section of the rail, above the sleeper, is shown in Figure 9 (b). Identically to the mid-span case, the rail response is quite similar for all viscoelastic models and the linear elastic model above the pinned-pinned frequency. Moreover, the first frequency peak for each viscoelastic model is similar to the mid-span case and the pinned-pinned anti-resonance frequency is unchanged.

Introduction of viscoelastic models in the ballast and in the pad

In this second study case, the pads are also modeled with a viscoelastic material. This comes from the fact that some rails have a configuration called slabtrack in which the rail and sleeper are linked by viscoelastic materials. Here, the Kelvin-Voigt model for the pad is used in the same time with various viscoelastic models for the ballast. The use of Kelvin-Voigt model is proposed for instance in [START_REF] Ntotsios | A comparison of ground vibration due to ballasted and slab tracks[END_REF]. The value of the damping constant which is calculated in these works and which is used here is 10 5 N.s.m -1 . The stiffness parameters for each direction are taken as the reference values from Table 4. Figure 10 (a) shows the vertical accelerance at mid-span. Introducing the Kelvin-Voigt model in the pad results in a variation of the vertical accelerance from the linear elastic model in the whole frequency range, notably above 1000 Hz. However, the pinned-pinned frequency remains unchanged. However, in a section above the sleeper, the pinned-pinned frequency changes and increases slightly as shows Figure 10 (b). The rail response above the pinned-pinned frequency is quite similar for all viscoelastic models. The influence of the complex moduli, as in the previous cases, is also observed at low frequencies.

Remarks on numerical results with the Duffing oscillator model in the ballast

Numerical problems solving (5) arise when a viscoelastic model is introduced in the pad. Singular Jacobian matrix and non-convergence occur when using the Kelvin-Voigt model or any other viscoelastic model reported in this work. In the case the pad is a linear elastic model, convergence is reached after 1746 iterations considering all frequency steps (CPU time about 0.3 s) for a criterion of 10 -6 in the relative error in the Newton-Raphson method. The initial solution a 0 for each frequency step is the converged solution of the last step.

General analysis : variation of accelerance due to viscoelastic damping in vertical direction

A general analysis of accelerance from the reference linear elastic model is eventually performed for each viscoelastic model at various longitudinal positions of a SAFE element. Let l = Lc be the longitudinal position in a SAFE element. In this analysis, c is taken from 0.5 to 1, this interval representing a position variation from mid-span to element extremity (above the sleeper). Due to the symmetry of the element, one does not need to perform calculations in the entire length span. In the following, the results are separated in three different application cases. The properties of the viscoelastic models are the same as those used in the previous sections.

Case 1: Application of viscoelastic models only in the ballast

In this first case, viscoelastic models are only applied to the ballast. The accelerance variation for the Zener and FD models are depicted in Figure 11 (a) and(b). Most of the accelerance increase appears at low frequency below 500 Hz and accelerance decrease is higher at specific frequencies such as 1000 Hz, 2000 Hz and 3000 Hz. The variation at those frequencies is uniform along the longitudinal position. Above 3000 Hz, the variation is close to zero. Similar variation patterns are observed for the GHM and GMM viscoelastic models and they are not represented here.

In order to study how accelerance changes in the overall SAFE element, the mean accelerance variation for each longitudinal position is computed over the frequency domain. Figure 12 (a) shows that it is almost constant with respect to c for each viscoelastic model. All viscoelastic models show an overall increase in accelerance, except the GHM model.

Case 2: Application of viscoelastic models in the ballast and Kelvin-Voigt model in the pad

In this second case, the viscoelastic models are applied to the ballast and a Kelvin-Voigt model for the pad. The accelerance variations for Zener and FD models are depicted in Figure 11 (c) and(d). Differently to the first case, accelerance increases mostly at low frequencies and increases at some specific frequencies such as 1000 Hz, 2000 Hz, 3000 Hz and 4500 Hz. However, these variations are no longer uniform along the SAFE element span. For instance, the variation is higher at 3000 Hz for positions in the interval 0.65 < c < 0.85 and at 4500 Hz in the interval 0.5 < c < 0.6 and 0.7 < c < 0.9. Similar results are observed for GHM and GMM models, they are not represented here.

Figure 12 (b) shows the mean accelerance variation with respect to c. An overall decrease in accelerance is observed for all viscoelastic models. Moreover, the accelerance reduction is more important at c = 0.95 and the slope is similar for all viscoelastic models. The GHM model brings more accelerance reduction than the other models.

Case 3: Application of viscoelastic models in the ballast and in the pad

In this third case, the viscoelastic models are applied to the ballast and to the pad in pairs (e.g Zener model at the ballast and the pad at the same time). The vertical accelerance variation in this case is shown in Figure 11 (e) and (f). Similarly to the first case, the accelerance increases mostly at low frequencies and is uniform along the SAFE element span. An important decrease is observed at 1000 Hz, 2000 Hz, 3000 Hz and 4500 Hz. The same behavior is observed for the GHM and GMM models. Figure 12 (c) shows the mean accelerance variation with respect to c. An overall increase in accelerance is observed, except for the GHM model.

Conclusions and prospects

In this work, we developed an accurate and cheap tool for analysing the effect of viscoelasticity in the rail vibration. We introduced linear and nonlinear viscoelastic behaviors in the ballast and in the pad in the SAFE method. The linear and nonlinear viscoelastic parameters were obtained either from the literature or by fitting the experimental data of storage modulus and loss modulus.

We analyzed the vertical accelerance of the top of the rail when a force is applied on it. When a viscoelastic model is applied only in the ballast, it is observed that the accelerance variation occurs mostly below the third frequency peak of the structure (1000 Hz). However, when a viscoelastic model is applied also in the pad, accelerance variation is observed in the whole frequency range.

We compared the SAFE method for some viscoelastic models applied in the ballast with a 1D Timoshenko beam model with periodic viscoelastic support. The comparison reveals that the beam model is not able to reproduce accelerances peaks because of the lack of cross-section deformation of the rail foot that appears in high frequency during a vertical excitation of the rail (foot flapping). In addition, in this work we proposed a solution for introducing the Duffing oscillator in the ballast. Using the harmonic balance method and the Newton-Raphson method, the numerical solution is successful when the Duffing oscillator is combined with linear elastic and Kelvin-Voigt pad. In this case, the results are satisfactory regarding the calculation time and numerical accuracy.

The main change in the code SAFE with respect to the linear elastic material is the implementation of the frequency-dependent complex stiffness of the viscoelastic materials, used for the calculation of the rail impedance. That operation is straightforward and no time-consuming. That model can be applied for the study and development of viscoelastic materials that can potentially mitigate the rolling noise.

As a perspective, we intend to implement some reduced order modeling (ROM) techniques. The ROM method is a powerful tool to provide rich information to designers, in terms of virtual charts that can be produced for various parameters of the system for further decision-making. It has already been applied in [START_REF] Barbarulo | Pgd-vtcr : A reduced order model technique to solve medium frequency broad band problems on complex acoustical systems[END_REF][START_REF] Barbarulo | Proper Generalized Decomposition applied to linear acoustic: a new tool for broad band calculation[END_REF][START_REF] Boucinha | Ideal minimal residual-based proper generalized decomposition for nonsymmetric multi-field models -Application to transient elastodynamics in space-time domain[END_REF][START_REF] Gonzalez | Real-time direct integration of reduced solid dynamics equations[END_REF][START_REF] Modesto | Proper generalized decomposition for parameterized Helmholtz problems in heterogeneous and unbounded domains: application to harbor agitation[END_REF][START_REF] Germoso | Efficient pgd-based dynamic calculation of nonlinear soil behavior[END_REF][START_REF] Borzacchiello | Advanced parametric space-frequency separated representations in structural dynamics : A harmonic-modal hybrid approach[END_REF][START_REF] Quaranta | From linear to nonlinear pgd-based parametric structural dynamics[END_REF] to elastodynamics problems. We intend to work with the PGD (Proper Generalized Method), in order to efficiently compute multi-parametrized solutions and conduct tractable optimization studies for the viscoelastic models.

Figure 2 .

 2 Figure 2. Standardized rail profile (here with norm 60E1), from [30] (left); Periodic 3D rail element, from[15] (right).

Figure 3 .

 3 Figure 3. Model for the rail and support components (pad, sleeper, ballast), from [11] (left); Zoom on the model of the periodic support with mass-spring system [15] (right).

  35.03 × 10 6 N/m, k N L = 4 × 10 14 N/m 3 and c L = 1732.50 × 10 3 Ns/m.

Figure 4 .

 4 Figure 4. (a) Storage modulus and (b) loss factor of the viscoelastic models.

Figure 5 .

 5 Figure 5. Rail geometry, mesh and force application and measure point.

Figure 6 .

 6 Figure 6. Vertical accelerance at mid-span and extremity (above the sleeper) of a SAFE element considering linear elastic constituents (the four frequency regions are depicted).

Figure 7 .

 7 Figure 7. Vertical accelerance response for a sleeper-viscoelastic ballast model

Figure 8 .

 8 Figure 8. Vertical accelerance of SAFE and beam model on periodic support at mid-span : (a) Zener, (b) FD, (c) GMM and (d) GHM.

Figure 9 .

 9 Figure 9. Viscoelastic model only in the ballast : Vertical accelerance at mid-span (a) and (b) Vertical accelerance above the sleeper.

Figure 10 .

 10 Figure 10. Viscoelastic model in the ballast and in the pad : (a) Vertical accelerance at mid-span and (b) Vertical accelerance above the sleeper.

Figure 11 .

 11 Figure 11. Vertical accelerance variation from the linear elastic case as function of span position c. Figures (a) and (b) are for Zener and FD models in case 1. Figures (c) and (d) are for Zener and FD models in case 2. Figures (e) and (f) are for Zener and FD models in case 3.

Figure 12 .

 12 Figure 12. Mean accelerance variation in the vertical direction : (a) Case 1, (b) Case 2 and (c) Case 3.

Table 1

 1 Linear elastic ballast and pad properties[START_REF] Baz | Active and passive vibration damping[END_REF]

Table 2

 2 Parameters of the GMM model

	Stiffness (N/m) Value Relaxation time (s) Value
	k 0	2.46x10 6	-	-
	k 1	0.247x10 9	ρ 1	2x10 -2
	k 2	0.247x10 9	ρ 2	2x10 -1
	k 3	0.347x10 9	ρ 3	2
	k 4	0.947x10 9	ρ 4	20
	k 5	0.347x10 9	ρ 5	200
	k 6	0.147x10 9	ρ 5	2000

Table 3

 3 Parameters of the GHM model with 3 oscillators.

	28.27 39559.1 8.99	4409	1.728	460.2	1	5.17x10 6

α 1 ω 1 (rad s -1 ) α 2 ω 2 (rad s -1 ) α 3 ω 3 (rad s -1 ) ζ (1,2,3) k 0 (N/m)

As a remark, we verified that the parameters used for the Duffing nonlinear model work in the nonlinear domain by setting k N L = 0 and we observed that the main effect is changing the accelerance at low frequency (below 400 Hz).