Benoît Orfeuvre
email: orfeuvre.benoit@gmail.com

Jérôme Tonetti

Gaël Kerschbaumer

Renaud Barthelemy

Alexandre Moreau-Gaudry

Mehdi Boudissa

Orfeuvre Benoît

EOS stereographic assessment of femoral shaft malunion after intramedullary nailing. A prospective series of 48 patients at 9 months' follow-up

Keywords: Femoral shaft fracture, Femoral nail, EOS, Malunion, Femoral version, Rotation disorder

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

Intramedullary nailing is the treatment of choice in femoral shaft fracture. One weakness is the risk of transverse, frontal or sagittal malunion. The most frequent deformity is transverse rotation disorder, with incidence of 20-30% [START_REF] Jaarsma | Rotational malalignment after fractures of the femur[END_REF][START_REF] Winquist | Closed intramedullary nailing of femoral fractures. A report of five hundred and twenty cases[END_REF][START_REF] Bråten | Femoral anteversion in normal adults. Ultrasound measurements in 50 men and 50 women[END_REF][START_REF] Yang | Prevention of malrotation deformity in femoral shaft fracture[END_REF].

To our knowledge, there have been few studies of malrotation in femoral shaft nailing, and measurement methods were poorly reproducible [START_REF] Jaarsma | Rotational malalignment after fractures of the femur[END_REF][START_REF] Winquist | Closed intramedullary nailing of femoral fractures. A report of five hundred and twenty cases[END_REF][START_REF] Bråten | Femoral anteversion in normal adults. Ultrasound measurements in 50 men and 50 women[END_REF][START_REF] Yang | Prevention of malrotation deformity in femoral shaft fracture[END_REF]. None used EOS stereoradiography [START_REF] Delin | Ionizing radiation doses during lower limb torsion and anteversion measurements by EOS stereoradiography and computed tomography[END_REF][START_REF] Kalifa | Evaluation of a new low-dose digital x-ray device: first dosimetric and clinical results in children[END_REF].

The aims of the present study were 1) to identify and quantify malrotation after transverse femoral shaft intramedullary nailing, using the EOS system; 2) to assess functional impact; and 3) to determine risk factors. The study hypothesis was that our surgical technique is associated with radiological rotation disorder rates equivalent to those reported in the literature.

Material and method

Population

A single-center continuous prospective study included patients aged >15 years, treated between January 2014 and December 2016 for unilateral femoral shaft fracture by anterograde intramedullary nailing. Joint fractures, floating knee and bilateral fractures were excluded, as were patients followed up elsewhere.

Ninety-six patients were eligible. Forty-eight were included: 12 female, (25%), 36 male (75%); mean age, 31.4 ± 14.8 years (range, 15.4-72.4 years). The inclusion flowchart is shown in figure 1. Mean body-mass index (BMI) was 23.6 ± 3.2 kg/m² (15.9-31,8). Fracture and surgery data are shown in Table 1.

Surgical technique

Patients were positioned on an orthopedic table with pubic support and boot traction of the foot. A titanium anatomic locking nail (T2, Stryker, Pusignan, France) was implanted.

Femoral neck anteversion was assessed in the transverse plane as the angle between the axis of the proximal femoral neck and the distal posterior bicondylar plane, as described by Bråten et al. [START_REF] Bråten | The role of fluoroscopy in avoiding rotational deformity of treated femoral shaft fractures: an anatomical and clinical study[END_REF]: fluoroscopy situated the femoral neck in relation to the ground after fracture reduction; strict lateral knee view superimposing the 2 condyles situated the bicondylar plane and thus the transverse axis of the knee, with respect to the ground; and subtracting the two axes gave the femoral neck anteversion value. The target value was 15-20° [START_REF] Bråten | The role of fluoroscopy in avoiding rotational deformity of treated femoral shaft fractures: an anatomical and clinical study[END_REF]. Varus-valgus (frontal plane) and genu recurvatum-flessum (sagittal plane) were assessed by aligning fragments according to the nailing principle. The entry point was at the summit of the greater trochanter. The reaming guide was positioned at the center of the condyles, frontally and laterally. Limb length restoration was checked by realigning the fragments.

After reaming and hardware introduction, static locking was performed proximally and freehand dynamic locking distally. Anteversion was adjusted, following Bråten, ahead of distal locking [START_REF] Bråten | The role of fluoroscopy in avoiding rotational deformity of treated femoral shaft fractures: an anatomical and clinical study[END_REF]. Immediate resumption of weight-bearing was at the surgeon's discretion.

Assessment

Patients were followed up at 45 days, 3 months, 6 months and 1 year, with standard AP and lateral femoral X-ray. At 6 months, clinical and radiological assessment including systematic EOS stereoradiography was performed by an independent observer. Data comprised epidemiological parameters, fracture characteristics, accident-to-surgery time, operator experience (senior or trainee), day or night surgery (emergency surgery after 6 pm), and clinical and radiological data comprising hip and knee ranges of motion, clinical measurement of morphological parameters including limb-length discrepancy, return to work, pain on VAS, and postoperative low back pain.

Radiologic assessment used EOS stereoradiography (EOS Imaging, Paris, France) [START_REF] Kalifa | Evaluation of a new low-dose digital x-ray device: first dosimetric and clinical results in children[END_REF][START_REF] Dubousset | EOS stereo-radiography system: whole-body simultaneous anteroposterior and lateral radiographs with very low radiation dose[END_REF]. The dedicated sterEOS software, version 1.6.4 (EOS Imaging, Paris, France) generated femur, tibia and fibula models for a 3D lower-limb skeleton model [START_REF] Dubousset | EOS stereo-radiography system: whole-body simultaneous anteroposterior and lateral radiographs with very low radiation dose[END_REF][START_REF] Chaibi | Fast 3D reconstruction of the lower limb using a parametric model and statistical inferences and clinical measurements calculation from biplanar X-rays[END_REF][START_REF] Than | Geometrical values of the normal and arthritic hip and knee detected with the EOS imaging system[END_REF]. The software automatically measured femoral torsion disorder: anteversion (°) with respect to the posterior bicondylar plane, limb-length discrepancy, valgus/varus, and flessum/genu recurvatum [START_REF] Chaibi | Fast 3D reconstruction of the lower limb using a parametric model and statistical inferences and clinical measurements calculation from biplanar X-rays[END_REF][START_REF] Than | Geometrical values of the normal and arthritic hip and knee detected with the EOS imaging system[END_REF].

On the anteroposterior axis, varus malrotation was defined by Hip-Knee-Ankle (HKA) angle <177° and valgus malrotation by HKA >183° [START_REF] Saragaglia | Computerassisted knee arthroplasty: comparison with a conventional procedure. Results of 50 cases in a prospective randomized study[END_REF][START_REF] Bellemans | The Chitranjan Ranawat Award: Is Neutral Mechanical Alignment Normal for All Patients?: The Concept of Constitutional Varus[END_REF]. On the mediolateral axis, normal alignment was defined by an 177-183° angle between the femoral axis (femoral head center to knee center) and the tibial axis (knee center to ankle center); angles >183° defined genu recurvatum compensating the malunion in flessum, and angles <177° defined flessum compensating malunion in genu recurvatum. There was assumed to be no anteroposterior translation, given the intramedullary nail design.

On the longitudinal axis, positive (≥15°) or negative (≤15°) differential axial rotation with respect to the healthy side defined malrotation (∆R), the main endpoint. Axial rotation was also analyzed as <0° or >0°, with respect to the surgical target of 15±5° (∆T). Limb-length discrepancy was measured as negative with respect to the contralateral length.

Stereoradiographic interpretation was validated by a single independent radiologist experienced in bone and joint imaging.

Statistics

Sample size for the main endpoint was 30 patients. Analyses used Stata 14. The study conformed to the Declaration of Helsinki, and had local IRB approval. All patients provided consent.

Results

Radiographic results

Fourteen patients (39.2%) showed malrotation (∆R): 11 with positive ∆R (mean, 27° ± 9.1 (range, 15-44)), 3 with negative ∆R (-23.1° ± 2.1 (-22 to -26°)) (figure 2).

The target (∆T) was not met for 32 patients (67%). Mean pathological side anteversion was 18.5° ± 13.8 (range, -17 to 63°). Twenty-four showed positive ∆T ( mean, 29.3° ± 8.0 (21-63)), and 8 negative ∆T (-2.6° ± 8.7 (-17 to 6) (figure 3).

Fifteen patients (31.2%) showed limb-length shortening (max., -13 mm) and 33 lengthening (68.8%) (max., +26 mm). Clinical limb-length discrepancy correlated with EOS measurements (rho 0.57; p<0.01).

On EOS, there was no significant difference between sides for frontal (p=0.22) or sagittal (p=0.42) rotation.

Clinical results

Clinical assessment was made at a mean 9.2 ± 5.1 months (range, 2.7-23.2) (Table 2).

Five patients (10.4%) resumed sports activity at 3 months, 12 (25%) at 3-6 months, 9 (18.8%) at 6-9 months, and 22 (45.9%) did not return to sport. None returned to competitive sport.

There was no clinical impact of limb-length discrepancy.

Complications

By end of follow-up, there had been 6 revision surgeries (12.5%): 1 correction of >20 mm limb-length discrepancy, 4 for non-union (at 7, 8, 11 and 13 months), and 1 femoral derotation procedure (at 4 months). There were 3 gluteus medius tendinites associated with bursitis. One patient showed symptomatic gluteus medius calcifications. Two showed signs of psoitis due to mechanical irritation by the proximal locking screw. There were no correlations between transversal, frontal or sagittal malrotation and SF12 or Harris scores (p>0.05). In the 11 patients with positive ∆R, mean SF12 was 7.6 ± 1.1 (range, 6-9), and mean Harris score 90 ± 8.9 (75-100). In the 3 patients with negative ∆R, mean SF12 was 7.3 ± 0.6 [START_REF] Bråten | The role of fluoroscopy in avoiding rotational deformity of treated femoral shaft fractures: an anatomical and clinical study[END_REF][START_REF] Dubousset | EOS stereo-radiography system: whole-body simultaneous anteroposterior and lateral radiographs with very low radiation dose[END_REF], and mean Harris score 71 ± 31 (36-96).

Multivariate analysis and radio-clinical correlations

Discussion

At a mean 9 months' follow-up, 14 of the 48 patients showed ≥15° malrotation on EOS with respect to the healthy side (29.2%). This is in agreement with the literature, with rates of 20-30% [START_REF] Jaarsma | Rotational malalignment after fractures of the femur[END_REF][START_REF] Fantry | Distal femoral complications following antegrade intramedullary nail placement[END_REF][START_REF] Hofmann | The role of intramedullary nailing in treatment of open fractures[END_REF].

One limitation of the present study was the large rate of loss to follow-up, due to patients changing residence. The use of EOS stereoradiography for assessment was also debatable, CT being the current gold-standard [START_REF] Jaarsma | Rotational malalignment after fractures of the femur[END_REF][START_REF] Jeanmart | Computer tomography of neck, chest, spine and limbs: atlas of pathologic computer tomography[END_REF]. The interest of EOS imaging to assess femoral torsion has been demonstrated [START_REF] Folinais | Measuring femoral and rotational alignment: EOS system versus computed tomography[END_REF][START_REF] Pomerantz | Three-dimensional biplanar radiography as a new means of accessing femoral version: a comparative study of EOS three-dimensional radiography versus computed tomography[END_REF], with no significant difference with respect to CT [START_REF] Pomerantz | Three-dimensional biplanar radiography as a new means of accessing femoral version: a comparative study of EOS three-dimensional radiography versus computed tomography[END_REF][START_REF] Escott | EOS low-dose radiography: a reliable and accurate upright assessment of lower-limb lengths[END_REF] or even superiority of EOS according to Morvan et al. [START_REF] Morvan | Femoral Torsion: Impact of Femur Position on CT and Stereoradiography Measurements[END_REF]. Radiation dose to the pelvis is 18.8fold lower than with classical X-ray [START_REF] Kalifa | Evaluation of a new low-dose digital x-ray device: first dosimetric and clinical results in children[END_REF], and, compared to CT, dose to the testicles is 24-fold lower, and 14-to-30-fold lower for the knee and ankle [START_REF] Delin | Ionizing radiation doses during lower limb torsion and anteversion measurements by EOS stereoradiography and computed tomography[END_REF]. The EOS system also enables simultaneous weight-bearing measurement of limb-length discrepancy, so that the healthy side can serve as reference. The limitations lie in interpretation and the evaluation system [START_REF] Knafo | Reproducibility of lowdose stereography measurements of femoral torsion after IM nailing of femoral shaft fractures and in intact femurs[END_REF]. For this reason, interpretation was made by an experienced independent radiologist. The evaluation software has been constantly upgraded since the report by Knafo et al. in 2016 [20].

Natural femoral anteversion is agreed to average 15°, with physiological variation between 0° and 30° [START_REF] Jaarsma | Rotational malalignment after fractures of the femur[END_REF][START_REF] Bråten | Femoral anteversion in normal adults. Ultrasound measurements in 50 men and 50 women[END_REF][START_REF] Bråten | Torsional deformity after intramedullary nailing of femoral shaft fractures. Measurement of anteversion angles in 110 patients[END_REF]. In the present study, mean postoperative anteversion was 18.5°, associated with a 29.2% rate of malrotation. This value was close to the 16° anteversion reported by Patel et al. after nailing for gunshot wound, but where the prevalence of malrotation was lower, at 12.3% [START_REF] Patel | Intramedullary nailing of diaphyseal femur fractures secondary to gunshot wounds: predictors of postoperative malrotation[END_REF]. In the present series, 67% of patients were off the target 15 ± 5°, suggesting that the method described by Bråten et al. is not very reproducible or reliable [START_REF] Delin | Ionizing radiation doses during lower limb torsion and anteversion measurements by EOS stereoradiography and computed tomography[END_REF].

There was a tendency for hyper-anteversion, with 50% of patients over the target, although without impact on SF12 or Harris functional scores.

There are several intraoperative methods intended to limit malrotation, without consensus [START_REF] Patel | Intramedullary nailing of diaphyseal femur fractures secondary to gunshot wounds: predictors of postoperative malrotation[END_REF][START_REF] Omar | Preoperative virtual reduction reduces femoral malrotation in the treatment of bilateral femoral shaft fractures[END_REF]. Some authors use preoperative external fixation, intraoperative CT control, or intraoperative navigation, with encouraging results [START_REF] Omar | Preoperative virtual reduction reduces femoral malrotation in the treatment of bilateral femoral shaft fractures[END_REF][START_REF] Ramme | Evaluation of malrotation following intramedullary nailing in a femoral shaft fracture model: Can a 3D c-arm improve accuracy?[END_REF][START_REF] Cao | Intraoperative Correction of Femoral Rotational Deformity Using a Conventional Navigation System and a Smartphone: A Novel Technique[END_REF].

There was a weak correlation between clinical assessment and EOS confirmation. All clinical methods are poorly reproducible, with up to 20° error in the hip [START_REF] Souza | Concurrent criterion-related validity and reliability of a clinical test to measure femoral anteversion[END_REF][START_REF] Nussbaumer | Validity and test-retest reliability of manual goniometers for measuring passive hip range of motion in femoroacetabular impingement patients[END_REF]. The present study found a correlation between internal rotation of the limb and femoral anteversion, an excess in one being associated with an excess in the other. Moreover, given a mean 30° physiological external tibial torsion, a phenomenon of adaptation to gait patterns in both pelvis and knee compensates malrotation [START_REF] Rosskopf | Femoral and tibial torsion measurement in children and adolescents: comparison of 3D models based on lowdose biplanar radiography and low-dose CT[END_REF]. Compensation for excess external rotation by the patient increases posterior acetabular and medial femorotibial stress, with risk of osteoarthritis of the hip and knee. Clinical findings of excess internal rotation indicate screening for malrotation following nailing. The absence of clinical impact of limb-length discrepancy exceeding 10 mm (max., +26mm) was unexpected [START_REF] Khamis | Relationship and significance of gait deviations associated with limb length discrepancy: A systematic review[END_REF], and may be due to short follow-up. Disturbed gait, partial functional recovery, recent trauma and hardware still in place are possible explanations; this requires longer-term analysis.

In the present study, the surgical revision rate was 12.5%, with 2.1% for femoral derotation.

The nonunion rate was 8.4% , in agreement with literature rates of up to 14.2% [START_REF] Malik | Factors affecting rates of infection and nonunion in intramedullary nailing[END_REF]. One reason could be soft-tissue incarceration in factures with large initial displacement, which is difficult to assess on standard X-ray or intraoperatively on the traction table. There were no cases of cortical breach or screw breakage.

  2 and R software version 3.4.1. The α risk was set at 5%. Normal distribution for quantitative variables was checked on Shapiro-Wilks test, and correlations were assessed on Pearson or Spearman correlation coefficient (rho) accordingly. Means were compared on Student test or nonparametric Mann-Whitney-Wilcoxon test. For qualitative variables with >2 modalities, ANOVA or non-parametric Kruskal-Wallis test was applied.

Table 3

 3 shows risk factors for malrotation: multi-site fracture (n=3) (p=0.04), surgeon inexperience (p=0.04), and open reduction (p=0.01). No factors emerged for limb-length discrepancy or frontal or sagittal malunion. Clinical and EOS femoral anteversion correlated (rho=0.53; p<0.01). Clinical internal rotation of the hip and EOS femoral anteversion correlated (rho=0.39; p<0.01).

  Risk factors for malrotation comprised: surgeon inexperience, or learning curve; multi-site fracture in multiple trauma; and open reduction. In our experience, open surgery does not prevent malrotation. A need for a sub-vastus approach on traction table testifies to difficulty in introducing the guidewire, and thus probably poor initial reduction. Associating traction and open reduction does not enable anatomic reduction following difficult intraoperative exposure.5. ConclusionEOSstereography assesses malrotation following closed nailing for femoral shaft fracture. Risk factors for malrotation comprised multi-site fracture, surgeon inexperience, and open reduction. Intraoperative bilateral comparative clinical and radiological assessment gives a better picture of native femoral version. This in-theater screening under fluoroscopy or CT can avoid certain severe rotation disorders in patients with natural hyperversion or retroversion. Disclosure of interest: The authors have no conflicts of interest to disclose in relation to the present study. Funding: none.

Figure 1

 1 Figure 1: Selection flowchart

Figure 2 :Figure 3 :

 23 Figure 2: Malrotation with respect to the healthy side (∆R). (Healthy side rotation = 0°)

Table 1 : Fracture characteristics and surgical data

 1

	Characteristics	Population, n (%) (N=48)
	Type of accident
	-road	28 (58.3)
	-sport	12 (25)
	-other	8 (16.7)
	Emergency surgery
	-no	6 (12.5)
	-yeas	42 (87.5)
	Surgeon's experience
	-junior	28 (58.3)
	-senior	20 (41.7)
	Open fracture
	-no	39 (81.2)
	-yes	9 (18.8)
	Side
	-right	17 (35.4)
	-left	31 (64.6)
	Type
	-unifocal	38 (79.2)
	-bifocal	7 (14.6)
	-multifocal	3 (6.2)
	Form
	-short spiroid	7 (14.6)
	-long spiroid	7 (14.6)
	-transverse	34 (70.8)
	Location on shaft
	-proximal	7 (14.6)
	-medial	39 (81.2)
	-distal	2 (4.2)
	Isolated osteoligamentous lesion
	-no	15 (31.2)
	-yes	33 (68.8)
	Multiple fracture
	-no	29 (60.4)
	-yes	19 (39.6)
	Open reduction
	-no	44 (91.7)
	-yes	4 (8.3)

Table 2 : Clinical results

 2

	Clinical parameters	N=48 (range)
	Femoral anteversion	10.6 ± 10.2° (0 -50)
	Scores
	-SF12 (/10 points)	7.2 ± 1.7 (3 -10)
	-Harris (/100 points)	85 ± 17.5 (36 -100)
	Pathologic hip range of motion
	-flexion	122.7 ± 15.9°
	-extension	-15.8 ± 8.9°
	-abduction	44.9 ± 6.2°
	-adduction	25.1 ± 11.1°
	-external rotation	50.9 ± 19.3°
	-internal rotation	32.0 ± 15.2°
	-pain on FABER	6 (12.5%)
	-pain on FADIR	8 (16.7%)
	Pathologic knee range of motion
	-flexion	143.6 ± 14.1°
	-extension	-1.1 ± 3.1°
	Meniscal pain
	-medial	6 (12.5%)
	-lateral	6 (12.5%)
	Knee extension
	-flessum (≥ 5°)	1 (2.1%)
	-genu recurvatum (-5 to -15°)	14 (29.2%)
	Frontal disorder
	-normal alignment	21 (43.8%)
	-varus	21 (43.8%)
	-valgus	6 (12.5%)
	Limb-length discrepancy
	-equality	20 (41.7%)
	-shortening	5 (10.4%)
	mean	-8 ± 6.7mm
	maximum	-20 mm
	-lengthening	23 (47.9%)
	mean	10.4 ± 4.9mm
	maximum	20 mm
	Acetabular pain syndrome	1 (2.1%)

Table 3 : Risk factors for intraoperative malrotation

 3

	Risk factors


Van-Vuong Lai
email: vanvuong.lai1990@gmail.com

Marc Anciant

Olivier Chiello

Jean-Franc ¸ois Brunel

Philippe Dufrénoy

A nonlinear FE model for wheel/rail curve squeal in the time-domain including acoustic predictions

Keywords: Curve squeal, Wheel/Rail contact, Finite element method, stability analysis, Time-domain (transient) analysis

Squeal noise of rail-bound vehicles frequently occurs in curves with a small radius and is a major nuisance for transport users and local residents. For the quantification of squeal intensity, a complete vibro-acoustic analysis is developed in this paper. This complete analysis requires time-domain analysis able to introduce non linearities leading to obtain dynamic saturation at the contact zone and a computation of sound radiation of the whole system. For time-domain analysis, a finite element (FE) formulation around the stationary position in an Eulerian reference frame is derived with a fine discretization of the contact surface combined with unilateral and Coulomb friction laws. Appropriate numerical techniques and reduction strategies are then used in order to solve the non linear discrete equations in dynamic self-sustained conditions. Both the transient approach and linear stability analysis are carried out. For sound radiation calculation, the contact forces calculated from wheel/rail contact model are then used for the calculation of squeal noise by using a coupled fluid-structure resolution based on boundary element method for the acoustic part and finite element method for the structural part. Results are first discussed in terms of unstable modes which are consistent between transient and stability analysis. Transient calculation shows that the apparent global friction coefficient during stick-slip cycles is slightly smaller than the constant local friction coefficient, and a dynamic saturation curve with hysteresis considerably different of the quasi-static curve. Finally the sound radiation calculation showed that the sound power radiated from the wheel is dominant with harmonics coming from the contact non linearities.

Introduction

Curve squeal is a type of railway noise which is highly occurred when tramway or subways come in a tight radius curve (radius lower than 200m [START_REF] Thompson | Railway noise and vibration: mechanisms, modelling and means of control[END_REF]). It is well accepted that the high lateral slip of the wheel on the rail-head is the main cause of curve squeal [START_REF] Rudd | Wheel/rail noise-part ii: Wheel squeal[END_REF][START_REF] Remington | Wheel/rail squeal and impact noise: What do we know? what don't we know? where do we go from here?[END_REF][START_REF] Vincent | Curve squeal of urban rolling stockpart 1: State of the art and field measurements[END_REF]. The friction forces generated by this sliding motion may lead to structural instability. Self-sustained vibration of the wheel/rail system is then developed leading to squeal. This noise is characterized by 5 high sound pressure levels (130 dB at 0.9 m from the wheel) at pure medium and high frequencies. The majority of measurements described by Rudd [START_REF] Rudd | Wheel/rail noise-part ii: Wheel squeal[END_REF], Vincent et al. [START_REF] Vincent | Curve squeal of urban rolling stockpart 1: State of the art and field measurements[END_REF], Koch et .al [START_REF] Koch | Curve squeal of urban rolling stockpart 2: Parametric study on a 1/4 scale test rig[END_REF] show that the wheel radiates at frequencies close to the natural frequencies of axial modes with zero nodal circle. However, squeal is occasionally observed at frequencies close to the natural frequencies of axial and radial wheel modes [START_REF] Glocker | Curve squealing of trains: Measurement, modelling and simulation[END_REF]. In addition, the highest noise level is often radiated by the inner leading wheel and rarely on the outer wheel [START_REF] Glocker | Curve squealing of trains: Measurement, modelling and simulation[END_REF]. Unlike rolling noise, the modeling of 10 curve squeal noise is less advanced because it has to take into account highly frequency dynamic events (stick/slip) at the contact zone during squeal. For a "quantitative" vibro-acoustic model, it requires firstly a time-domain analysis, which can describe nonlinear events leading to obtain dynamic saturation at the frictional rolling contact zone and then a computation of sound radiated from the wheel and the rail.

Contact models, which allow for the description of phenomena at wheel/rail contact zone, are essential for all curve 15 squeal models. It exists several models based on analytic formulas or heuristic laws allowing for the determination of contact forces. In quasi-static condition (no dynamic), they are Hertz theory [START_REF] Hertz | Über die berührung fester elastischer körper[END_REF] for the normal problem and Kalker linear theory [START_REF] Kalker | Wheel-rail rolling contact theory[END_REF], Vermeulen and Johnson model [START_REF] Vermeulen | Contact of nonspherical elastic bodies transmitting tangential forces[END_REF], Shen et al. model [START_REF] Shen | A comparison of alternative creep force models for rail vehicle dynamic analysis[END_REF] or Chopaya law proposed by Ayasse et al. [START_REF] Ayasse | 4 wheel-rail contact[END_REF] for the tangential problem. In dynamic condition, Gross-Thebing [START_REF] Gross-Thebing | Lineare Modellierung des instationären Rollkontaktes von Rad und Schiene[END_REF] proposed a law in which the low frequency creep force -creepage law is represented by a damper to which a spring is introduced in series to simulate 20 the high frequency transient effects. For rolling noise, Hertz theory and Gross-Thebing's law are used and perform well in comparison with experimental results [START_REF] Thompson | Experimental validation of the twins prediction program for rolling noise, part 1: description of the model and method[END_REF][START_REF] Thompson | Experimental validation of the twins prediction program for rolling noise, part 2: results[END_REF]. However, these laws are only applied for relatively small amplitudes and creepage. These quasi-static models may not be sufficient well to estimate curve squeal characterized by high nonlinear phenomena (stick/slip) and high creepage occurring at the contact zone. On the other hand, by using Kalker's discretized contact models such as simplified theory [START_REF] Kalker | A fast algorithm for the simplified theory of rolling contact[END_REF] in [START_REF] Xie | Introduction of falling friction coefficients into curving calculations for studying curve squeal noise[END_REF][START_REF] Huang | Squeal prediction for a bogied vehicle in a curve[END_REF][START_REF] Ding | Effects of rail dynamics and friction characteristics on curve squeal[END_REF] or variational theory [START_REF] Kalker | Wheel-rail rolling contact theory[END_REF] in [START_REF] Pieringer | A numerical investigation of curve squeal in the case of constant wheel/rail friction[END_REF][START_REF] Zenzerovic | An engineering time-domain model for curve squeal: Tangential point-contact model and green's functions approach[END_REF][START_REF] Zenzerovic | Influence of spin creepage and contact angle on curve squeal: A numerical approach[END_REF][START_REF] Giner-Navarro | Study of railway curve squeal in the time domain using a high-frequency vehicle/track interaction model[END_REF], the impact of some simplifications (Boussinesq and Cerruti elastic half-space assumptions, contact/friction decoupling) rather unknown in the case of high frequency dynamics and very fast evolutions at the contact zone. The modeling of wheel/rail rolling contact by the finite element (FE) method has been recently developed. The majority of wheel/rail FE contact models are performed within a Lagrangian frame [START_REF] Wiest | Assessment of methods for calculating contact pressure in wheel-rail/switch contact[END_REF][START_REF] Zhao | The solution of frictional wheel-rail rolling contact with a 3d transient finite element model: Validation and error analysis[END_REF][START_REF] Toumi | Finite element analysis of the frictional wheel-rail rolling contact using explicit and implicit methods[END_REF][START_REF] Yang | Numerical modeling of wheel-rail squeal-exciting contact[END_REF]. They are in good agreement with Hertz theory and CONTACT software based on Kalker's variational theory [START_REF] Kalker | Wheel-rail rolling contact theory[END_REF] in both normal and tangential solutions in a quasistatic state. However, nonlinear transient FE simulations with frictional contact still remain very expensive in terms of CPU time and memory size because the potential contact area in the rolling direction needs to be meshed with elements of sufficient small size to insure the desired precision. To reduce the computational time, a finite element (FE) computational method for the dynamics of frictional rolling contact systems in an Eulerian reference frame was developed in [START_REF] Lai | Full finite element models and reduction strategies for the simulation of friction-induced vibrations of rolling contact systems[END_REF]. Continuous equations of the problem are derived around the stationary position of rolling in an Eulerian reference frame. This approach allows for the large model reduction because the potential contact zone become very small in comparison with those in Lagrangian frame. Moreover, in order to obtain reasonable computation times, reduction strategies were also proposed in [START_REF] Lai | Full finite element models and reduction strategies for the simulation of friction-induced vibrations of rolling contact systems[END_REF]. This model finally allows both transient and stability analyses to predict unstable frequencies and full spectrum of friction-induced vibrations. The results of the quasi-static analysis show a good agreement with the ones obtained with Kalker's variational theory (implemented in CONTACT software) and analytical laws. Concerning the performance of the reduction strategies, the contact static approximation performs well for the reduction of computational time but still ensures the desired precision. In an unstable configuration, results of the time-domain analysis provided solutions which are coherent with the stability results. Unlike contact models or friction characteristic, less attention has been paid for the computation of wheel/rail sound radiation in curve. Considering that the sound radiated from the rail is very smaller than the sound radiated from the wheel in the squeal phenomenon [START_REF] Vincent | Curve squeal of urban rolling stockpart 1: State of the art and field measurements[END_REF], most of curve squeal models only calculate the noise radiated from the wheel [START_REF] Huang | Squeal prediction for a bogied vehicle in a curve[END_REF][START_REF] Zenzerovic | An engineering time-domain model for curve squeal: Tangential point-contact model and green's functions approach[END_REF][START_REF] Zenzerovic | Influence of spin creepage and contact angle on curve squeal: A numerical approach[END_REF][START_REF] Schneider | Noise generation in railway wheels due to rail-wheel contact forces[END_REF][START_REF] Fingberg | A model of wheel-rail squealing noise[END_REF][START_REF] Périard | Wheel-rail noise generation: curve squealing by trams[END_REF]. Schneider and Popp [START_REF] Schneider | Noise generation in railway wheels due to rail-wheel contact forces[END_REF] evaluated the radiated sound using the Rayleigh integral for a baffled plate. Fingberg [START_REF] Fingberg | A model of wheel-rail squealing noise[END_REF], Huang [START_REF] Huang | Squeal prediction for a bogied vehicle in a curve[END_REF] and Zenzerovic [START_REF] Zenzerovic | An engineering time-domain model for curve squeal: Tangential point-contact model and green's functions approach[END_REF][START_REF] Zenzerovic | Influence of spin creepage and contact angle on curve squeal: A numerical approach[END_REF] used a combination of the boundary element method (BEM) and modal expansion techniques. By using this technique, the sound power at a frequency is evaluated as a sum of contributions from wheel modes. Techniques using in these models allow for the estimation of the sound radiated by the wheel. However, without the computation of the sound radiated by the rail, the contribution of rail to this noise has been still unknown.

The aim of this paper is to develop a FE model for quantifying curve squeal. Considering that curve squeal is a consequence of self-sustained friction-induced vibrations of the wheel-rail system, in the first part, a nonlinear timedomain analysis allowing for the determination of these vibrations is developed in comparison with linear stability analysis. These analyses are based on the finite element contact model developed in [START_REF] Lai | Full finite element models and reduction strategies for the simulation of friction-induced vibrations of rolling contact systems[END_REF]. The most performant reduction basis including free-interface modes and residual static attachment modes (CSA) is used to reduce the computational times. Models such as [START_REF] Rudd | Wheel/rail noise-part ii: Wheel squeal[END_REF] consider that friction coefficient decreasing as a function of creep velocities can generate instabilities of wheel/rail contact, acting as a negative damping introduced by the slope of creep/friction curve. However recent works have shown that instabilities could also be obtained with a constant friction coefficient, associated to mode coupling [START_REF] Glocker | Curve squealing of trains: Measurement, modelling and simulation[END_REF][START_REF] Pieringer | A numerical investigation of curve squeal in the case of constant wheel/rail friction[END_REF] and damping behavior of the track [START_REF] Ding | Effect of rail dynamics on curve squeal under constant friction conditions[END_REF][START_REF] Lai | The critical effect of rail vertical phase response in railway curve squeal generation[END_REF]. Thus, it is introduced constant at the local scale. Results are discussed and compared with those obtained with the point-contact model for both stability and time-domain analysis. In the second part, the computation of sound radiation for both wheel and rail is performed using the results of the time-domain analysis and based on the boundary element method (BEM). The contribution of wheel/rail to this noise is evaluated quantitatively.

Models for wheel/rail rolling contact and acoustic radiation

Statement of the problem

A wheel/rail rolling contact is described in Fig. 1. The z-axis coincides with the common normal to the two surfaces in contact, the longitudinal x-axis corresponds to the rolling direction and the y-axis refers to the lateral direction. The wheel rolling speed is V.

In curve, three movements of wheel may happen [START_REF] Kalker | Wheel-rail rolling contact theory[END_REF][START_REF] Ayasse | Determination of the wheel rail contact patch in semi-hertzian conditions[END_REF]: lateral sliding, longitudinal sliding or rotation around their common normal (or spin) with relative velocities ∆V y , ∆V x and ∆Ω z respectively. ∆V y can be approximated by: Vα where α is angle of attack. As a result of these movements, the relative sliding instantaneous velocity (or creep velocities) between the wheel and the rail at the contact interface in the Eulerian frame is given by [START_REF] Johnson | Contact mechanics[END_REF] : where subscripts R, W denote respectively the rail and the wheel. The notation u = ∂u ∂t refers to the time partial derivative. u(x, y, z, t) and v(x, y, z, t) denote respectively the displacement and velocity fields of the structure. The terms involving V ∂u ∂x represent the deformation contributions due to rolling in the Eulerian frame whereas the terms involving u simply represent the dynamic contributions.

ṡx = v W x -v R x = ∆V x -∆Ω z y + V( ∂u W x ∂x - ∂u R x ∂x ) + (u W x -uR x ) ṡy = v W y -v R y = ∆V y + ∆Ω z x + V( ∂u W y ∂x - ∂u R y ∂x ) + (u W y -uR y ) (1)

Transient and stability analysis

In order to evaluate the response of wheel/rail contact system excited by these creep velocities (Eq. ( 1)), the finite element contact model developed and validated in [START_REF] Lai | Full finite element models and reduction strategies for the simulation of friction-induced vibrations of rolling contact systems[END_REF] is adopted. Both transient and stability analysis are performed. Local constant friction is assumed at the contact interface. For the time-domain analysis, non-smooth laws such as Signorini's unilateral law and Coulomb's law to deal with frictional contact on the interface. These laws are written in equivalent semi-regularized forms ( [START_REF] Jean | The non-smooth contact dynamics method[END_REF]). Convective terms are neglected. The principle of virtual power is used for the system dynamics with contact laws written in weak forms. The integration scheme called "modified θ-method" [START_REF] Loyer | Study of nonlinear behaviors and modal reductions for friction destabilized systems. application to an elastic layer[END_REF] which allows to introduce nonlinearities (impact, stick, slip) is used to compute transient solutions. An iterative fixed point algorithm on equivalent contact reactions and friction forces is used for the non linear resolutions for quasi-static and dynamic solutions at each time step.

It must be noted for the following equations that M, K, C are the mass, stiffness and viscous damping matrix of the whole system. U, R n and R t denote respectively the vectors of nodal displacements, normal reactions and friction forces. These nodal displacements, normal reactions and friction forces depend on the relative sliding instantaneous velocity in Eq. ( 1) and friction coefficient. Their determination is presented in [START_REF] Lai | Full finite element models and reduction strategies for the simulation of friction-induced vibrations of rolling contact systems[END_REF]. In order to reduce the size of the system and the computational time for the nonlinear resolutions at each time step, the "Contact static approximation" which performs well for the reduction of computational time but still ensures the desired precision [START_REF] Lai | Full finite element models and reduction strategies for the simulation of friction-induced vibrations of rolling contact systems[END_REF] is chosen in this study. The principle is to search an approximated solution U = Bq r of the problem spanned by a reduced basis B, leading to a reduced dynamics equation. B is a composed of free-interface normal modes Φ enriched by static residual attachment modes Φs which can be written:

Φs = Φ s -Φ(Φ T KΦ) -1 Φ T P c T ( 2
)
where P c T = [P n T P t T ] and P n , P n are matrices allowing to pass the contact reactions from the local relative frame to the global frame. Φ s is the static displacement responses to unit contact reactions e.g. attachments modes such as:

KΦ s = P c
T . Hence, Φs corresponds to attachment modes after the elimination of the contribution of normal modes. Using basis B = [Φ Φs ] together with the elimination of the dynamic terms relating to the residual attachment modes Φs , this reduction strategy consists of two steps: global step by solving the global dynamics using free-interface normal modes (Eq. ( 3)) and local step by adding a local static residual flexibility controlled by matrix ΦT s (Eq. ( 4)) in the expression of the contact displacements:

Φ T MΦ q + Φ T CΦ q + Φ T KΦq = Φ T (F + P n T R n + P t T R t ) (3) ΦT s K Φs q s = ΦT s (F + P n T R n + P t T R t ) (4)
where q and q s are the generalized coordinate vectors corresponding respectively to normal modes and residual attachment modes such that U = Φq + Φs q s . In parallel with time-domain analysis, stability analysis is carried out around full sliding equilibriums in order to predict unstable modes. The linearization of the non-linear equations around the equilibrium leads to a constrained non symmetric eigenvalues problem. Complex modes are then computed. Modes corresponding to complex eigenvalues λ with positive real part are unstable. The divergence rate of a complex mode is defined as (λ)/ (λ) where (λ), (λ) are respectively the real and imaginary parts (pulsation) of the mode. The same reduction technique is used for the stability analysis but considering residual attachment modes Φs defined on the effective contact region at the quasi-static equilibrium, leading to the following eigenvalue problem in Eq. ( 5):

λ 2 Φ T MΦ + λΦ T CΦ + Φ T (K + K c )Φ q = 0 (5)
where

K c = ( PT n + µ PT t ) I n ΦT s ( PT n + µ PT t ) -1
Pn is a non-symmetrical matrix taking into account the effects of the frictional forces. µ is friction coefficient. Pn , Pt are new projections matrices on nodes in the effective contact region at equilibrium. I n is the boolean localization matrix such that

Pn = I n Pc with Pc = [ PT n PT t ]
. This stability analysis then consists in solving the global dynamics using free-interface modes in Eq. ( 5) coupled by a local static residual stiffness controlled by matrix K c .

Calculation of the squeal noise radiated by a railway system: wheel/rail/sleepers

The methodology used for the calculation of the sound radiated by both wheel and track has been developed by ESI-Group. Compared to rolling noise simulation [START_REF] Thompson | Railway noise and vibration: mechanisms, modelling and means of control[END_REF], where the contact forces are related to the roughness of the rail and the wheel contact surfaces, in the case of squeal noise, the origin of the forces is the alternation of (stick/slip) phases in the contact zone, which can be predicted by a nonlinear mechanical simulation in the time domain as described before.

The vibro-acoustic simulation uses these contact forces as an input to compute, in the frequency domain, the acoustic power radiated by one of the three system: (1) the wheel, (2) the two rails and (3) the sleepers.

Assuming the three systems to be linear, the acoustic power radiated by one system (i), noted W i a , can be expressed as a quadratic form of the two contact forces ( f y, f z), such as:

W i a = [ f y f z ] * H i ayy H i ayz H i azy H i azz [ f y f z ] T (i=1,2,3...: wheel, rails, sleeper...) (6)
where H i a is the (2x2) matrix containing acoustic transfer functions, with unit (W/N 2 ) and () * states for complex conjugate operation. In this study, the extra-diagonal terms of the acoustic transfer matrix that represent the crosscorrelation terms between acoustic pressure and velocity are neglected. This means that the power radiated is the summation of the power due to the lateral force with that due to the vertical component. Also, it is assumed that the vibration on the three systems wheel/rails/sleepers are uncorrelated one from the other, and therefore the total acoustic power W a is given by the simple summation: W a = (i=1:3) W i a . Force spectrum is obtained by using Fourier transformation for stationary temporal forces (limit cycle) as presented in section 3.4.3 For each system, the acoustic power is computed by using a coupled fluid-structure approach based on a Boundary Element Method (BEM) [START_REF] Hamdi | Development of a sound radiation model for a finite-length duct of arbitrary shape[END_REF] for the acoustic part and a Finite Element Method (FEM) for the structural part. Considering the wheel for example, it is meshed with solid elements for the structural analysis, whereas the BEM mesh is a surface mesh corresponding to the external skin of the wheel. The degrees of freedom in the simulation are respectively the acoustic pressure (P) at the nodes of the BEM mesh and the displacement vector w = (w x , w y , w z ) at the nodes of the structural mesh, and these are solution of the following symmetric linear system [START_REF] Hamdi | A mixed functional for the numerical resolution of fluid-structure interaction problems[END_REF]

K -ω 2 M -ρω 2 A B + C/2 B T + C T /2 -D/(ρω 2 ) w P = F -F a S a (7)
where K, M are the stiffness and mass matrices of the structure. D(ω) is the acoustic admittance BEM matrix. A(ω) is the pseudo-added mass matrix due to the fluid. B(ω) is the BEM coupling matrix. C is the FEM coupling matrix. S a , F a : acoustic source and force vectors. F : mechanical force vector. ω = 2π f is the pulsation of the problem with f the frequency (Hz).

A modal approach also is used for the structure in order to reduce the size of the system. Once the linear system in Eq. ( 7) is inversed, the acoustic power is obtained by performing the following surface integral on the BEM mesh:

W i a = S 1/2 (PV n )dS ( 8
)
where P is the acoustic pressure field and V n the normal velocity field, both defined on the external skin of the wheel. The (normal) velocity field on the skin of the structure is classically related to the (normal) displacement by:

V n = iωw n . (•) denotes the real part of (•). This vibro-acoustic simulation also gives the opportunity to extract (for the excitation at hand) the radiation efficiency indicator, which tells about the capacity of the system to transform vibration into acoustic power. This indicator, noted σ, which usually tends to the unity value for high frequencies, is given by:

σ = W a S ρC V n 2 with V n 2 = 1 2S S | V n 2 |dS (9)
3. Application to wheel/rail dynamics in curve

In this section, an application of the numerical method presented in the previous sections is proposed for a wheel/rail model in rolling contact with imposed lateral creepages as in a curve squeal situation.

Wheel/Rail models

The wheel/rail models considered in [START_REF] Lai | The critical effect of rail vertical phase response in railway curve squeal generation[END_REF] are used in this paper. The wheel FE model corresponds to a standard steel monobloc wheel with nominal rolling diameter of 920 mm and a mass of 314 kg (Fig. 2). From this model, 100 free-interface modes have been calculated up to 8000 Hz considering clamped boundary conditions at the hub. 3 kinds of mode may be distinguished: the radial modes (r, n), the axial modes (a, m, n) and the circumferential modes (c, m) where n is the number of nodal diameters and m is the number of nodal circles [START_REF] Thompson | Railway noise and vibration: mechanisms, modelling and means of control[END_REF]. Modal damping factors are chosen depending on the nodal diameters (1, 0.1 and 0.01% for n = 0, n = 1 and n ≥ 2 respectively, cf. [START_REF] Thompson | Railway noise and vibration: mechanisms, modelling and means of control[END_REF]). The track FE model consists in one periodically supported rail of UIC60 type as used in [START_REF] Lai | The critical effect of rail vertical phase response in railway curve squeal generation[END_REF]. This rail is 48 m long. The space between the sleepers (sleeper span) is 60 cm. The FE track model is made of 200000 quadratic elements and 600000 DOFs. As a first approach, the dynamics of the sleepers and the ballast is neglected because of its rather low frequency domain. The rail has 48 m length but ends with 2 anechoic terminations of 6 m composed of 5 rail portions of length L = (0.6; 0.6; 0.6; 1.2; 3.0) m with increasing damping, avoiding the return of waves. The track support contains only elastic pads that connect the rail and each sleeper. They are modeled by 69 springs of longitudinal, lateral and vertical stiffnesses (K x = Ky = 36, Kz = 180) MN/m for each sleeper. The contact position is in the center of the rail in the x-axis. The rail structural damping is η = 0.02. An equivalent viscous damping model is derived from structural damping factors of 0.02 for the rail, 1 for the pad and 0.1 to 1 for the rail anechoic terminations [START_REF] Betgen | On the prediction of rail cross mobility and track decay rates using finite element models[END_REF]. For both wheel/rail models, the contact zone is meshed with elements size of 1mm (Figs. 3 and4). To validate this mesh, result is compared with that obtained from analytical models in quasi-static condition [START_REF] Lai | Full finite element models and reduction strategies for the simulation of friction-induced vibrations of rolling contact systems[END_REF].

Results of the stability analysis

In this section, the results of the stability analysis of the wheel/rail system for the reference case is firstly presented. The bifurcation curves (effect of friction coefficient) then allow to analyze the influence of the discretization of the contact zone on instabilities.

Reference case analysis

The reference case with the kinematic parameters listed in Table 1 is considered in this section. Results of the stability analysis is carried out in case of full lateral sliding using a higher creepage s yo = ∆V y /V = 1%. In order to solve the non symmetric eigenvalue problem, the reduction strategy presented in Eq. ( 5) (contact static approximation) is performed. Four unstable modes with positive real part are obtained as listed in Tab. 2. The corresponding unstable mode shapes are shown in Fig. 5. As for the point-contact model [START_REF] Lai | The critical effect of rail vertical phase response in railway curve squeal generation[END_REF], there are 3 axial wheel modes with (2,3,4) nodal diameters and zero nodal circle involved in the first three unstable modes and some axial and radial wheel modes involved in the fourth unstable mode. The predictions are in good agreement with general experimental results [START_REF] Vincent | Curve squeal of urban rolling stockpart 1: State of the art and field measurements[END_REF] which showed that the main squeal frequencies are close to natural frequencies of axial modes with zero nodal diameter.

The corresponding frequencies are almost the same between the two models whereas the divergence rates are higher for the surface-contact model.

Bifurcation curves

The bifurcation curves of the unstable complex modes are represented in Fig. 6 in comparison with the results obtained with the point-contact model [START_REF] Lai | The critical effect of rail vertical phase response in railway curve squeal generation[END_REF]. The critical friction coefficients above which the complex modes become unstable are respectively equal to 0.06, 0.1, 0.22 and 0.3. Fig. 6 shows also that with or without discretization of the contact area, the unstable frequencies are almost identical (difference < 2Hz). However, the divergence rates of the unstable modes are higher for the surface contact model than for the point contact model due to the difference between the Hertz's stiffness k H in [START_REF] Hertz | Über die berührung fester elastischer körper[END_REF] and the FE local stiffness matrix (K c in Eq. ( 5). An "equivalent" average FE stiffness k f e can be calculated, equal to the sum of all terms of (I n ΦT s ( PT n + µ PT t )) -1 . Stiffness k f e is found to be 1.19 times larger than Hertz's stiffness k H . On the other hand, the form of the bifurcation curves of each unstable mode obtained with the surface-contact model is similar to that obtained with the point-contact model. It can be conclude that the instability mechanisms are similar for the two models. As found in [START_REF] Lai | The critical effect of rail vertical phase response in railway curve squeal generation[END_REF], the first three unstable modes are due to equivalent damper behavior of the rail whereas the fourth unstable mode results from a mode coupling instability. The above results show that for stability analysis, the discretization of the contact zone does not modify substantially the instability mechanisms. However due to a stronger coupling between the normal contact DOFs, the divergence rates of the unstable modes are greater for the surface-contact model than for the point-contact model.

Results of the time-domain analysis

From the stability analysis of the wheel/rail contact performed in the previous section, four unstable modes are obtained. However, stability analysis does not allow to introduce nonlinearities and determine the amplitude and full spectrum of the friction-induced vibration. The aim of this section is to perform time-domain analysis in order to compare these analyses, especially in the case of several unstable modes and to test the influence of some parameters.

Reference case analysis

Results of the time-domain analysis corresponding to the unstable case founded in the previous section are determined using a numerical time integration from given initial conditions. In all the following results, the integration starts from the equilibrium i.e. the initial displacements are the displacements obtained from the quasi-static solution and the initial velocities are null. The parameters V = 10 m/s, µ = 0.3 and s yo = ∆V y /V = 1% are used. The time step for the integration is ∆t = 1µs.

The time series of the lateral velocity of a point outside the contact zone is presented in Fig. 7a. The spectrogram in Fig. 7b allows to observe which frequencies are present in the solution and when they appear. At the beginning of the simulation, there are two major frequencies: 334 and 919 Hz which are very close to the natural frequencies of the two first unstable modes provided by stability analysis. In the stationary step, where nonlinear events happen, the dominant frequency is f o = 919 Hz and its harmonic frequencies f k = k f o appear. This dominant frequency corresponds to the unstable mode with the greatest positive real part (Tab. 2). Fig. 8a shows the time series of the lateral contact resultant forces F t . The tangential resultant force increases until a pronounced creep/slip oscillation builds up as shown in Fig. 8b. When the tangential contact resultant force is smaller than the traction bound µF n , a transient stick zone appears at the leading edge of the effective contact region as shown in Fig. 9. The apparent global friction coefficient which is defined as a ratio between the resulting tangential and normal contact force is represented in Fig. 10. This coefficient varies between 0.2958 corresponding to a partial slip state and 0.3 corresponding to the full slip state. It is clear that the average global friction coefficient during these stick-slip cycles is lightly smaller than the constant local friction coefficient. Thus, the decrease of the global friction coefficient could be interpreted as a consequence (and not a cause) of the friction-induced vibrations. 1) representing the deformation contributions due to rolling in the Eulerian frame, the global dynamic lateral creepage can be defined as:

0 0.
s yd = s yo + 1 V ∂u y ∂t ( 10
)
where ∂u y ∂t is the averaged lateral dynamic velocity in the effective contact zone. The dynamic saturation curve is defined as the relation between the dynamic tangential resulting contact force and the global dynamic lateral creepage. The curve is represented in Fig. 11 in comparison with some quasi-static point-contact models. An hysteresis curve is obtained with the full FE model. This curve is considerably different of quasi-static saturation curve. A possible explanation is that the quasi-static models do not take in account the lateral contact stiffness which is a key parameter at high frequencies [START_REF] Thompson | Railway noise and vibration: mechanisms, modelling and means of control[END_REF][START_REF] Knothe | Modelling of railway track and vehicle/track interaction at high frequencies[END_REF][START_REF] Mindlin | Compliance of elastic bodies in contact[END_REF]. These results of the this reference case show a consistence with the results of the stability analysis and exhibit localized stick/slip oscillations in the contact zone as in [START_REF] Pieringer | A numerical investigation of curve squeal in the case of constant wheel/rail friction[END_REF]. The influence of parameters on final limit cycles are analyzed in following sections.

Influence of imposed lateral creepage

Time-domain analysis is performed with s yo % varying from -0.9% to 0.9% while other parameters in Tab. 1 as well as initial conditions are maintained constant. Friction-induced vibrations are only obtained with s yo >= 0.5% (Tab. 3), which corresponds to the full sliding state of the contact zone between inner wheel and rail. Besides, results of the stability analysis does not show instabilities for contact between outer wheel and rail (s yo < 0) . These results are in good agreement with experimental observations in [START_REF] Vincent | Curve squeal of urban rolling stockpart 1: State of the art and field measurements[END_REF] which showed that squeal is mainly observed near the inner wheel. s yo -0.9% -0.8% -0.7% -0.5% -0.4% 0.4% 0.5% 0.7% 0.8% 0.9% friction-induced vibrations 0 0 0 0 0 0 x x x x Table 3: results of the time-domain analysis obtained with difference imposed lateral creepage (0: no friction-induced vibrations, x: stick-slip limit cycle)

Influence of friction coefficient

To study the influence of friction coefficient, time-domain analysis is performed with friction coefficient µ = 0.4 and other parameters in Tab. 1 as well as initial conditions are maintained constant. For reference purposes, stability analysis gives the unstable frequencies and corresponding type of unstable modes shapes in Tab. 4.

UF (Hz) Real part

Wheel modes 334 7.3 (a,0,2) 917.9

13.9 (a,0,3) 1670 2.05 (a,0,4) 3417

Table 4: Stability results obtained with surface-contact model for µ = 0.4 (UF: unstable frequency)

The time series and the spectrogram of the lateral velocity of a point outside the contact zone is presented in Fig. 12. At the beginning of this simulation, there are three major frequencies whose amplitude increases: 330, 920 and 3430 Hz which are very close to the frequencies of the three unstable modes with the greatest real part. In the second step where nonlinear events happen, there is competition between these modes. The first mode (330Hz) stops its contribution in the solution. The two other modes continue to develop. In the final step, the frequency f o = 3430 Hz and its harmonic frequencies f k = 2 f o become dominant. This dominant frequency corresponds to the unstable mode with the greatest positive real part as shown in Tab. 4. From the results in these above sections, it seems that the dominant frequency in the limit cycle corresponds to the unstable mode with the greatest real part. However, Loyer et al [START_REF] Loyer | Study of nonlinear behaviors and modal reductions for friction destabilized systems. application to an elastic layer[END_REF] showed that in the case with several unstable modes, different limit cycles can be found depending on initial conditions. On the other hand, for the modeling of curve squeal, little existing models have discussed about the influence of initial conditions on final limit cycle. This influence has been considered (not detailed in the paper) and exhibits that, also for the modeling of curve squeal, initial conditions act on the limit cycle and dominant frequencies. It also modifies the amplitude of the saturation curve but does not change the profile of the dynamic response. In terms of geometry, the FE model used for the vibro-acoustic simulation is similar to the one used previously for the mechanical squeal simulation, but the level of detail requested is not as high especially in the contact area. The average mesh size (∼5cm for the track, ∼1.5cm for the wheel) is ruled by a classic quarter wavelength criterion. It enables to properly model the wave propagation up to 6kHz for both systems. The FE model of the track presented Fig. 13 is made of 2 rails and 46 sleepers. It is 21.6m long. Due to the symmetry of the track along the rail X-axis, only half of the track is modeled and special boundary conditions are applied in order to properly simulate the complete track behavior, whose initial length is 43.2m. The track is attached to the ground by a set of damped springs and a 6 meters anechoic termination has been added in order to minimize the reflecting waves in the rails. The mesh contains 245988 nodes and 55560 solid elements. The acoustic mesh corresponds to the skin of the track. It is obtained by a surface coating operation and it contains 42666 nodes and 81882 surface elements. An infinite rigid acoustic plane is added bellow the track in order to simulate the wave reflection of the rigid ground made of concrete. When exciting each system individually with a unit force in vertical or lateral direction, this enables to extract the acoustic transfer function, defined as the ratio between the acoustic power radiated by the system, and the squared force amplitude as defined in Eq. ( 8). Fig. 14a depicts the acoustic transfer functions obtained for the track (rails and sleepers) considering the vertical force only, and Fig. 14b the acoustic transfer functions of the wheel for both load cases lateral and vertical. For this excitation, the power radiated by the rails is clearly higher than the one radiated by the sleepers, once the decoupling frequency ( f ≈ 300 Hz) between rails and sleepers is passed. This cut-off frequency is driven by the stiffness of the springs placed between the rails and the sleepers and the masses at hand. On Fig. 15a is compared the radiation efficiency of the rail obtained with the current BEM-coupled 3D method, and the analytical result obtained for an equivalent infinite cylinder with a radius R= 8.6 cm moving as a rigid body in lateral direction [START_REF] Thompson | Railway noise and vibration: mechanisms, modelling and means of control[END_REF]. A good correlation between both results is observed especially at low frequencies.

For the wheel (Fig. 14b), the radiated power in both directions exhibits lots of sharp peaks at eigenfrequencies of the wheel, whose amplitude is only controlled by the structural damping. The damping values retained for the simulation are quite weak (10 -3 , 10 -4 ) and depends on the azimuth order of the mode, as suggested in reference [START_REF] Thompson | Railway noise and vibration: mechanisms, modelling and means of control[END_REF]. Bellow 1kHz, the lateral force clearly produces more acoustic power than the vertical one. The radiation efficiency of the wheel depicted Fig. 15b, shows values close to one around the first axial mode (a,0,0) near 228 Hz, but also for frequencies greater than 500 Hz. It can be noticed that the mode at 336 Hz does not radiate that much, because for circumferential modes with even order (m=2,4,..) acoustic contributions tend to cancel each other. The results in terms of transfer function (in particular the peaks corresponding to resonances for the wheel) and of radiation efficiency (radiation efficiency close to 1 from 1000 Hz) are not surprising. The same results can be found as for rolling noise [START_REF] Thompson | Railway noise and vibration: mechanisms, modelling and means of control[END_REF] since they are independent of the excitation.

Squeal acoustic power results

The nonlinear mechanical analysis for the reference case performed previously produces squeal contact forces as time signals ( f y (t), f z (t)). These have to be transformed from the time domain into the frequency domain in order to be combined with the acoustic transfer functions previously computed. This is done by performing a Fourier transform of the time signals, with a proper normalization, resulting in force Power Spectral Densities (PSD) ( f 2 y (ω), f 2 z (ω)), with unit (N 2 /Hz). Fig. 16a represents an example of force time signals ( f y (t), f z (t)) for a flexible rail configuration with imposed lateral creepage s yo = 1%. The graph represents a time slot of T = 0.01 s, but the total signal duration processed is T = 0.1 s and the sample frequency is F s = 10 6 Hz. Both signals obviously look periodic. The force level in vertical direction is 13dB higher than the one in lateral direction, but less harmonic distortion is observed for this signal which looks like a pure tone. In the frequency domain, the force PSDs presented Fig. 16b self-sustained vibrations due to the unstable dynamic behavior of mechanical systems with friction such as brake squeal [START_REF] Sinou | A global strategy based on experiments and simulations for squeal prediction on industrial railway brakes[END_REF][START_REF] Soobbarayen | Noise and vibration for a self-excited mechanical system with friction[END_REF][START_REF] Corradi | Prediction of squeal noise based on multiresolution signal decomposition and wavelet representationapplication to fem brake systems subjected to friction-induced vibration[END_REF], the excitation due to curve squeal is clearly different from that of the rolling noise [START_REF] Thompson | Railway noise and vibration: mechanisms, modelling and means of control[END_REF]: there is a fundamental frequency which corresponds to a resonance and harmonics which do not correspond to resonances of the wheel. The fundamental frequency is close to unstable frequency found in stability analysis, which is also close to natural frequency of the wheel. These harmonic frequencies result from nonlinear events during squeal. The frequency range of the excitation being mainly greater than 1000 Hz, one may almost use a radiation efficiency equal to 1 and do without the vibro-acoustic calculation to make an approximate estimate of the radiated power levels. The simple multiplication of the force PSD by the acoustic transfer functions (and the summation for both forces) as reported in Eq. ( 7) leads to the acoustic power PSD (W/Hz) for the different systems. The strong excitation level at fundamental frequency 920Hz and harmonics are retrieved as peaks in the final power PSD spectra presented in Fig. 17a. The global acoustic power level L=121.3dB is fully dominated by the wheel. At 920Hz, the high response level can be obviously explained by the fact the force excitation level at that frequency is high, but also by the fact the order of the wheel mode excited (m=3) is odd, and therefore the transformation of the vibration into acoustic power by the wheel is even more efficient. At 1837 Hz, no such coincidence occurs and the second peak in the response simply corresponds to the second harmonic in the excitation signal. For a lower imposed lateral creepage s yo = 0.5%, the force excitation levels are 15-20 dB lower: (Fy, Fz) = (19.5, 36.8) dBN for s yo = 0.5% whereas (F y , F z ) = (39.6, 52.5) dBN for s yo = 1%, but the general shape of the force PSD are still the same, with a fundamental at 920 Hz plus harmonics. The acoustic power radiated by the wheel (Fig. 17b) consequently is also ( 17 dB) lower: L=103.9dB against L=121.3 dB.

Conclusion

In this paper, a model for wheel/rail frictional rolling contact in curve is developed. For rolling contact model, a FE model developed and validated in [START_REF] Lai | Full finite element models and reduction strategies for the simulation of friction-induced vibrations of rolling contact systems[END_REF] is adopted. The reduction basis including the free-interface modes and static residual flexibility modes is used in order to reduce the computation loads and times. Time-domain analysis is also performed in order to determine the amplitude and full spectrum of the self-sustained vibrations in curve.

results of the time-domain analysis are used to compute the sound radiation of wheel and rail. Stability analysis is performed in comparison with the results obtained with the point-contact model in order to highlight the influence of the discretization of the contact zone on the instability mechanism and intensity.

The stability results show that the discretization of the contact zone does not modify the instability mechanisms. However, due to a stronger coupling between the normal contact DOFs, the divergence rates of the unstable modes are greater for the surface-contact model than for the point-contact model.

The results of the transient calculation are consistent with the results of the stability analysis and exhibit localized stick/slip oscillations in the contact zone. The apparent global friction coefficient during stick-slip cycles is slightly smaller than the constant local friction coefficient. Thus, the decrease of the apparent global friction coefficient could be considered as a consequence of the friction-induced vibrations. The dynamic saturation curve with hysteresis obtained by the full FE model is considerably different of the quasi-static curve. This is probably due to the lateral wheel/rail contact stiffness at high frequencies. Thus, it is necessary for the curve squeal point-contact models to take into account this lateral stiffness. Even when these predictions are in good agreement with general experimental results in terms of squealing wheel or modes, the data will have to be validated experimentally.

Although the goal of the rail dynamic play an essential role for instabilities of wheel/rail contact system as found in recent works [START_REF] Ding | Effects of rail dynamics and friction characteristics on curve squeal[END_REF][START_REF] Lai | The critical effect of rail vertical phase response in railway curve squeal generation[END_REF][START_REF] Ding | An assessment of mode-coupling and falling-friction mechanisms in railway curve squeal through a simplified approach[END_REF], the sound radiation calculation shows that the sound power radiated from the wheel is clearly dominant especially at the peak corresponding to the fundamental frequency. Thus, the computation of the sound radiated by the wheel is enough sufficient to evaluated the squeal sound radiation.

Figure 1 :

 1 Figure 1: Coordinate system, kinetic variables (α: angle of attack of the curve, ∆V y : imposed lateral sliding velocity, ∆V x : imposed longitudinal sliding velocity and ∆Ω z : imposed spin velocity)

Figure 2 :

 2 Figure2: Wheel FEM mesh (wheel of type "Vyksa BA005" with a nominal rolling diameter of 920 mm and a mass of 314 kg) with a rigid constraint applied at the inner face of the hub, where the wheel is connected to the axle.

Figure 3 :

 3 Figure 3: FE model of the wheel with fine mesh on the contact zone

Figure 4 :

 4 Figure 4: FE model of the rail with fine mesh on the contact zone

( a )Figure 5 :

 a5 Figure 5: Unstable mode shapes

Figure 6 :

 6 Figure 6: Bifurcation curves of the unstable modes obtained with the FE model (solid line) and point-contact model (dotted line)

Figure 7 :Figure 8 : 9 -

 789 Figure 7: Tangential velocity of a point outside the contact zone

Figure 9 :

 9 Figure 9: Status of nodes in the potential contact zone: no contact zone (red), slip zone (blue) and stick zone (green)

Figure 11 :

 11 Figure 11: Dynamic saturation curve. The contact status at points marked with Arabic numerals is represented Fig. 9

Figure 12 :

 12 Figure 12: Tangential velocity of a point outside the contact zone for µ = 0.4

3. 4 .

 4 Sound radiated by a railway set wheel/rails/sleepers 3.4.1. Numerical FE models for vibro-acoustic simulation

Figure 13 :

 13 Figure 13: Finite Element model of the track used for vibro-acoustic simulation

Frequency

  exhibit a fundamental frequency at f = 920 Hz and regularly spaced harmonics at frequencies: 1840, 2760, 3680 Hz,... Like

Figure 14 :Frequency

 14 Figure 14: Acoustic transfer functions (W/N 2 ) for the track (rails and sleepers) and the wheel in the frequency range [80 -5630] Hz, a) for vertical force and b) for vertical and lateral forces

Figure 15 :

 15 Figure 15: Radiation efficiency of the rail (a) obtained with BEM simulation and equivalent cylinder analytic model extracted from [1] and of the wheel (b) obtained with BEM-coupled simulation for lateral and vertical loads

ForceFigure 16 :

 16 Figure 16: Force time signals f y (t) and f z (t) as well as their PSDs, ( s yo = 1%)

FrequencyFigure 17 :

 17 Figure 17: Acoustic power PSD for wheel, rails and sleepers

Table 1 :

 1 Kinematic parameters of the wheel/rail rolling contact model

			Rolling velocity V (m/s)		10
		Imposed lateral velocity of the wheel V y (m/s)	0.1
			Friction coefficient µ		0.3
			Static vertical load N	70000 N
		Point-contact model		Surface-contact model
	UF (Hz) Real part	Wheel modes	UF (Hz) Real part	Wheel modes
	334	4.0	(a,0,2)	334	5.11	(a,0,2)
	918.3	8.0	(a,0,3)	917.9	9.55	(a,0,3)
	1671	1.0	(a,0,4)	1670	1.0	(a,0,4)
	3418	1.38

(a,0,6), (r,5) and (a,1,1) 3417 2.73 (a,0,6), (r,5) and (a,1,1)

Table 2 :

 2 Stability

results obtained with both point-contact and surface-contact models in the reference case (UF: unstable frequency)

Acknowledgement

This work was carried out within the CERVIFER project backed by the ADEME organization and the Haut-de-France region and within the framework of the LABEX CeLyA (ANR-10-LABX-0060) of Université de Lyon, within the program " Investissements d'Avenir " (ANR-11-IDEX-0007) operated by the French National Research Agency (ANR). The authors gratefully acknowledge the support of these institutions.