Valérie Demange
email: valerie.demange@inrs.fr.

Michel Grzebyk

Michel Hery

Nicole Massin

Christophe Paris

Michel Héry

Pascal Wild

Longitudinal predictors of bronchial hyperresponsiveness and FEV1 decline in bakers

Keywords:

à la diffusion de documents scientifiques 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.

A c c e p t e d m a n u s c r i p t INTRODUCTION Baker's asthma is one of the most frequent occupational asthma (OA) in industrialized countries [START_REF] Jeebhay | Food processing and occupational respiratory allergy-An EAACI position paper[END_REF][START_REF] Page | Exposure to flour dust and sensitization among bakery employees[END_REF]). The prevalence is estimated to be from 4 to 13 % [START_REF] Jeebhay | Food processing and occupational respiratory allergy-An EAACI position paper[END_REF]. In France, flour was still the most frequent causative agent of OA during the period 2008-2014 [START_REF] Iwatsubo | Observatoire National des Asthmes Professionnels II. Bilan de la phase pilote et perspectives pour la surveillance des asthmes en lien avec le travail[END_REF]). Baker's asthma has seen the list of its causative agents expand over the years with the evolution of baking techniques, with the use of cereal flour, fungal enzymes, other ingredients like white egg proteins but also the presence of storage mites, arthropods, and molds [START_REF] Quirce | Diagnosis and management of grain-induced asthma[END_REF].

Bronchial hyperresponsiveness (BHR) is considered the physiological hallmark of asthma, whether occupational or not [START_REF] Malo | An official American Thoracic Society Workshop Report: presentations and discussion of the fifth Jack Pepys Workshop on Asthma in the Workplace. Comparisons between asthma in the workplace and non-work-related asthma[END_REF] [START_REF] Lipinska-Ojrzanowska | Bronchial Response to High and Low Molecular Weight Occupational Inhalant Allergens[END_REF]). Thus, it is a common hypothesis that BHR is a stage on the pathway to asthma [START_REF] Malo | An official American Thoracic Society Workshop Report: presentations and discussion of the fifth Jack Pepys Workshop on Asthma in the Workplace. Comparisons between asthma in the workplace and non-work-related asthma[END_REF][START_REF] Radon | Career Advice for Young Allergy Patients[END_REF]. A longitudinal follow up of the BHR status in populations exposed to asthmogens seems therefore to be relevant.

Apart in asthma, BHR is also observed in other lung diseases, such as chronic obstru ctive pulmonary disease [START_REF] Borak | Bronchial hyperresponsiveness[END_REF]. BHR was shown to be a predictor of irreversible loss of lung function in asthmatics with or without pulmonary symptoms [START_REF] Borak | Bronchial hyperresponsiveness[END_REF], but also in asymptomatic individuals from a Swiss population-based cohort [START_REF] Brutsche | Bronchial hyperresponsiveness and the development of asthma and COPD in asymptomatic individuals: SAPALDIA cohort study[END_REF]) and in followed-up firefighters [START_REF] Aldrich | Bronchial Reactivity and Lung Function After[END_REF]). However, BHR was not predictive of lung function loss in a recent relatively short -term study among alumin iu m smelter workers [START_REF] Abramson | Predictive value of non-specific bronchial challenge testing for respiratory symptoms and lung function in aluminium smelter workers[END_REF]).

An association of occupational flour dust exposure with a decreased lung function has been found in flour mill [START_REF] Bagheri | Relationship between lung function and flour dust in flour factory workers[END_REF][START_REF] Demeke | Assessment of Respiratory Symptoms and Pulmonary Function Status among Workers of Flour Mills in Addis Ababa, Ethiopia: Comparative Cross-Sectional Study Pulm Med[END_REF][START_REF] Lagiso | Chronic respiratory symptoms, lung function and associated factors among flour mill factory workers in Hawassa city, A c c e p t e d m a n u s c r i p t southern Ethiopia: "comparative cross-sectional study[END_REF][START_REF] Mohammadien | Effects of exposure to flour dust on respiratory symptoms and pulmonary function of mill workers Egypt[END_REF][START_REF] Zamani | The effect of chronic exposure to flour dust on pulmonary functions[END_REF] or bakery workers [START_REF] Fahim | Pulmonary function impairment and airway allergy among workers in traditional bakeries[END_REF][START_REF] Moghaddasi | Health-risk assessment of workers exposed to flour dust: a cross-sectional study of random samples of bakeries workers[END_REF][START_REF] Said | Effects on respiratory system due to exposure to wheat flour Egypt[END_REF], as recently reviewed by Zamani [START_REF] Zamani | The effect of chronic exposure to flour dust on pulmonary functions[END_REF]. In these studies, the occupational flour dust exposure was solely represented by the exposure status [START_REF] Fahim | Pulmonary function impairment and airway allergy among workers in traditional bakeries[END_REF][START_REF] Lagiso | Chronic respiratory symptoms, lung function and associated factors among flour mill factory workers in Hawassa city, A c c e p t e d m a n u s c r i p t southern Ethiopia: "comparative cross-sectional study[END_REF][START_REF] Moghaddasi | Health-risk assessment of workers exposed to flour dust: a cross-sectional study of random samples of bakeries workers[END_REF][START_REF] Mohammadien | Effects of exposure to flour dust on respiratory symptoms and pulmonary function of mill workers Egypt[END_REF], the duration of exposure [START_REF] Demeke | Assessment of Respiratory Symptoms and Pulmonary Function Status among Workers of Flour Mills in Addis Ababa, Ethiopia: Comparative Cross-Sectional Study Pulm Med[END_REF][START_REF] Said | Effects on respiratory system due to exposure to wheat flour Egypt[END_REF], or by the atmospheric flour dust concentration at the time of the data collection [START_REF] Bagheri | Relationship between lung function and flour dust in flour factory workers[END_REF][START_REF] Zamani | The effect of chronic exposure to flour dust on pulmonary functions[END_REF]. Different spirometric parameters were affected depending on the study [START_REF] Zamani | The effect of chronic exposure to flour dust on pulmonary functions[END_REF]. All these studies were cross sectional, with the exception of the two-year follow-up study by Zamani et al. Indeed, many of the published longitudinal studies in bakers concerned rather apprentices [START_REF] Acouetey | Genetic predictors of inflammation in the risk of occupational asthma in young apprentices Annals of allergy, asthma & immunology : official publication of the American College of Allergy[END_REF][START_REF] De Meer | Bronchial responsiveness to adenosine-5'-monophosphate and methacholine as predictors for nasal symptoms due to newly introduced allerge ns. A follow-up study among laboratory animal workers and bakery apprentices[END_REF][START_REF] De Zotti | Prospective study of work related respiratory symptoms in trainee bakers[END_REF][START_REF] Demange | Do airway inflammation and airway responsiveness markers at the start of apprenticeship predict their evolution during initial training? A longitudinal study among apprentice bakers, pastry makers and hairdressers[END_REF][START_REF] Gautrin | Incidence and determinants of IgE-mediated sensitization in apprentices. A prospective study[END_REF][START_REF] Skjold | The incidence of respiratory symptoms and sensitisation in baker apprentices[END_REF][START_REF] Tossa | Increase in exhaled nitric oxide is associated with bronchial hyperresponsiveness among apprentices[END_REF][START_REF] Walusiak | The risk factors of occupational hypersensitivity in apprentice bakers --the predictive value of atopy markers[END_REF]. These studies A c c e p t e d m a n u s c r i p t brought knowledge about the time course of the work-related respiratory symptoms or diseases at the beginning of exposure. The highest incidence of these symptoms seems to be in the first 6 to 12 months of employmen t [START_REF] Radon | Career Advice for Young Allergy Patients[END_REF]. The recognised risk factors are atopy [START_REF] Skjold | The incidence of respiratory symptoms and sensitisation in baker apprentices[END_REF], pre-existing sensitization to occupational allergens, [START_REF] De Meer | Bronchial responsiveness to adenosine-5'-monophosphate and methacholine as predictors for nasal symptoms due to newly introduced allerge ns. A follow-up study among laboratory animal workers and bakery apprentices[END_REF]BHR (de Meer et al. 2003;[START_REF] Demange | Do airway inflammation and airway responsiveness markers at the start of apprenticeship predict their evolution during initial training? A longitudinal study among apprentice bakers, pastry makers and hairdressers[END_REF][START_REF] Radon | Career Advice for Young Allergy Patients[END_REF][START_REF] Skjold | The incidence of respiratory symptoms and sensitisation in baker apprentices[END_REF]). We identified four longitudinal studies in bakers at work [START_REF] Cullinan | Allergen and dust exposure as determinants of work-related symptoms and sensitization in a cohort of flour-exposed workers; a case-control analysis[END_REF][START_REF] Florentin | Exhaled nitric oxide and screening for occupational asthma in two at-risk sectors: bakery and hairdressing[END_REF][START_REF] Kim | The usefulness of job relocation and serum eosinophil cationic protein in baker's asthma[END_REF][START_REF] Remen | Early incidence of occupational asthma is not accelerated by atopy in the bakery/pastry and hairdressing sectors[END_REF]. None of these studies were focused on BHR and lung function evolution, except for Kim et al., who focused on BHR changes after 6 month job relocation in a small population of 12 bakers.

The present longitudinal study aimed at describing the natural history of BHR and FEV1 decline by identifying their respective predictors, in a population of bakers and a group of non-exposed controls, the latter providing a baseline for the exploration of the dose-response relationships

MATERIALS AND METHODS

Study design and inclusion criteria

The study design was a longitudinal follow-up of the workforces at 4 industrial bakeries, with visits in 1992, 1994, 1996, 2001, 2003 and 2005 (n=110) and of a non-exposed group of blue-collar workers (n=38) at a salt-packing factory with visits in 1993, 2002, 2004 and 2006. Each visit comprised a clinical examination, a standardized respiratory questionnaire, occupational allergen skin prick tests (SPT), pulmonary function tests and a methacholine bronchial challenge test, all carried out at the workplace. The cross-sectional analysis of the initial data collected in 1992 was previously published [START_REF] Bohadana | Respiratory symptoms and airway responsiveness in apparently healthy workers exposed to flour dust[END_REF]. The baker jobs were croissant making , dough mixing, special breads making, general bread making, control and packaging, deliveries, cleaning, oven handling. The non-exposed salt-packing factory workers were employed in production, maintenance, or quality control departments. All workers who volunteered to participate were included. However, workers with an acute respiratory infection were not included in the visit. The non-exposed subjects were included if they had not been occupationally exposed to flour dusts or to another respiratory health hazard. As the aim of this paper was to determine long-term predictors, we included in the analyses all subjects who had attended at least two visits. Figure 1 presents the flowchart of the participation of the subjects in the study.

Respiratory health status

Respiratory symptoms were recorded at each medical examination by the same experienced interviewer using the European Coal and Steel Community questionnaire [START_REF] Bohadana | Respiratory symptoms and airway responsiveness in apparently healthy workers exposed to flour dust[END_REF][START_REF]Questionnaire of the European Community for Coal and Steel (ECSC) on respiratory symptoms. 1987--updating of the 1962 and 1967 questionnaires for studyi ng chronic bronchitis and emphysema[END_REF]). The questionnaire included detailed smoking habits, history of respiratory diseases and childhood allergic disea ses such as asthma and eczema. Clinical atopy was defined as having childhood allergic diseases (asthma or eczema). At each visit, we abstracted respiratory symptoms defined as a wheezing chest or a runny nose in absence of any cold.

Spirometry was carried out by the same 3 experienced technicians using the same electronic spirometer (Spiro -Analyzer ST 300, Fukuda Sangyo Co. Tokyo, Japan) in compliance with the European Respiratory Society recommendations [START_REF] Quanjer | Lung volumes and forced ventilatory flows. Report Working Party Standardization of Lung Function Tests, European Community for Steel and Coal. Official Statement of the European Respiratory[END_REF]. The observed / predicted FEV1 ratio was used for description, using the reference equations from Quanjer [START_REF] Quanjer | Lung volumes and forced ventilatory flows. Report Working Party Standardization of Lung Function Tests, European Community for Steel and Coal. Official Statement of the European Respiratory[END_REF].

We used 3 additive doses of methacholine delivered by a nebulizer (Mediprom FDC 88 -Paris, France), according to an abbreviated version of the methacholine bronchial challenge (MBC) test [START_REF] Gardner | Report of Snowbird workshop on standardization of spirometry[END_REF]. This corresponds to a maximal inhalated dose of 11.5 µmol, i.e. 2 100 µg. The first dose of 0.5 µmol is equivalent to one puff fro m an aerosolized solution with a concentration of 0.36 mg/ml, using a nebulizer delivering 0.14 ml/minute during 2 minutes. The second dose of 2.5 µmol is equivalent to one puff from an aerosolized solution with a concentration of 1.79 mg/ml. The third and final dose of 7.5 µmol is equivalent to one puff from an aerosolized solution with a concentration of 5.36 mg/ml.

The challenge test was discontinued either after the inhalation of the last dose of methacholine or if the FEV1 fell by 20% or more below the baseline value obtained prior to MBC testing. In that case, subjects were classified as having a positive MBC test, i.e. having bronchial hyperresponsiveness (denoted BHR in the sequel).

Skin prick tests to occupational allergens

Three SPTs were performed using wheat flour, rye flour and alpha amylase respectively (Stallergenes Laboratories, Fresnes, France). A 9% codeine phosphate solution was used as a positive control and the 50% glycerinated diluent as a negative control. A positive SPT was defined as a wheal diameter at least 3 mm larger than those obtained with the negative controls after 20 minutes. Sensitization to occupational allergens was defined at each visit as a positive reaction to at least one of the three allergens tested.

Exposure assessment

In each bakery, some full-shift dust samples of the inhalable fraction were obtained in order to assess the exposure of each job assignment. Closed-face filter holders were used (Millipore, Massachusetts, USA), housing pre-A c c e p t e d m a n u s c r i p t weighed 37-mm diameter polyvinyl chloride filters (Pall Gelman Sciences, Champs sur Marne), connected to portable battery-operated vacuum pumps (Escort, MSA, Châtillon sur Chalaronne) sampling at airflow rates of 1 l/mn. The PVC filters were weighed twice: in the laboratory during their preparation before sampling, and again in the same laboratory after sampling. These weighings were done in a controlled and constant atmosphere for heat and humidity, after a stay of at least twelve hours in the room. Before the first weighing, the filters were treated with Triton to eliminate electrostatic charges. The pumps were calibrated before and after sampling.

Generally, the samples were taken on workers, with the filter holder fixed near the breathing zone. In exceptional cases stationary sampling was used: for team leaders in charge of administrative tasks in a room isolated from the workshop. They were reserved for workers present in this room in a static position and who did not intervene directly in the manufacturing process. The pumps and sample holders were placed a few tens of centimeters away from their upper respiratory tracts. This measurement of the airborne exposure could thus be considered as representative of their individual exposure.

For a given job, the exposure level was calculated as the mean value of the dust concentration measurements obtained for this job. As these measurements could not be considered an exact measurement of the chronic exposure, these exposure means were categorized into three groups, <1mg/m3, 1-5 mg/m3, >5mg/m3 which were considered reliable by the experienced industrial hygienist (MH) who was part of the research group. This exposure level was assigned to all workers of the job. Exposure measurements were repeated several times over the entire study period, at least each time the working conditions (new equipment, automation of a productio n line, etc.) changed.

Statistics

According to the main objective of this study, we studied the respective determinants of BHR and respiratory symptoms, and the determinants of the evolution of FEV1 during the study follow-up.

Firstly, we studied whether BHR, respiratory symptoms and their joint occurrence were associated with gender, cumulative smoking, clinical atopy, sensitization to occupational allergens and exposure to flour dust (level and duration) using a mixed logistic regression with the subject as a random factor to account for the longitudinal nature of the data. Secondly, the repeated FEV1 measurements were modelled using linear mixed models with the subject as random effect. Several models were considered, including cumulative smoking, pre-visit and concurrent BHR, sensitization to occupational allergens , duration and exposure level, as independent variables . All analyses The analysis was carried out with STATA statistical software (STATA, College Station, TX, USA).

RESULTS

One hundred and ten bakers and 38 non-exposed subjects attended the inclusion visit and at least one follow-up visit, representing 483 visits. That is 3.3 visits per subject with a mean duration of follow-up of 6 years among the bakers and 8 years among the non-exposed subjects. Table 1 gives a basic description of subjects at inclusion, from inclusion to last follow-up, and at the last visit. At inclusion, clinical atopy was less often reported by bakers (11%) than non-exposed workers (21%), whereas sensitization to occupational allergens was more often observed in bakers (20% versus 16%). At the last visit, this prevalence was 30% in bakers versus 10% in non -exposed workers. At inclusion, the exposed workers had been working as bakers for 10.6 years on average and the most frequent category of exposure level to inhalable dust was between 1 to 5 mg/m 3 at inclusion and at the last visit. 18.8 (15.4) 16.5 (14.3) a : defined as the average of the inhalable dust concentrations for the different jobs held by the worker Table 2 describes the health outcomes of interest according to the exposure status, at inclusion and the last visit.

At inclusion, the prevalence of BHR was 24% in bakers and increased to reach 36% at the last visit. In non -exposed workers, this prevalence increased slightly between these two visits (8 and 11%, respectively). Moreover, in bakers at both visits, the prevalence of BHR increased with duration of exposure. Such a pattern was not observed with the exposure level. The prevalence of respiratory symptoms was similar at inclusion in exposed and non -exposed workers (38% versus 34%, respectively), while it was higher in bakers at the last visit (52% versus 40%, respectively). For this outcome, no pattern with duration of exposure or mean level of exposure was observed. The observed FEV1 was close to the predicted value in bakers both at inclusion and the last visit, and slightly higher in non-exposed workers. In bakers, at both visits, the FEV1 in percent predicted decreased with duration of exposure.

As for BHR prevalence, no pattern was observed with the exposure level. : defined as fall in FEV1 of 20% or more during the MBC test; b : defined as a wheezing chest or a runny nose in absence of any cold; c : observed / predicted FEV1 ratio using the reference equations from [START_REF] Quanjer | Lung volumes and forced ventilatory flows. Report Working Party Standardization of Lung Function Tests, European Community for Steel and Coal. Official Statement of the European Respiratory[END_REF].

A c c e p t e d m a n u s c r i p t

The probability of having BHR (Table 3) doubled when the duration of exposure increased by ten years (adjusted OR=2.2, 95% CI=[1.0;5.0]). Female gender and the number of pack-years for ever-smokers also significantly increased the probability of having BHR. The probability of reporting respiratory symptoms was multiplied by almost 4 when the duration of exposure increased by ten years (adjusted OR=3.6, 95% CI=[1.4;9.4]). Clinical atopy was significantly associated with respiratory symptoms (adjusted OR=16.5, 95% CI=[1.7;157.1]). When considering having BHR and at least one respiratory symptom simultaneously, duration of exposure and clinical atopy were significant predictors. The exposure level had no statistical significant effect on any of the three outcomes considered, neither did a history of sensitization to occupational allergens. Table 4 presents the results of two multiple linear regression analyses of ΔFEV1, the difference between observed FEV1 and its predicted value [START_REF] Quanjer | Lung volumes and forced ventilatory flows. Report Working Party Standardization of Lung Function Tests, European Community for Steel and Coal. Official Statement of the European Respiratory[END_REF] according to gender, height and age. In the first model, ΔFEV1

decreased by almost 200 ml per ten years of exposure ). BHR at a preceding visit was significantly (although borderline) associated with ΔFEV1 even after adjustment for the significant effect of BHR at the current visit. Similarly to the results for symptoms and BHR, there was no effect of the exposure level. Finally, as expected, ΔFEV1 decreased significantly with smoking (by 6 ml per pack year).

The second model excluding the (non-significant) exposure level and including the history of occupational sensitization confirmed the relationship with duration of exposure, the association with BHR at the current visit and cumulative smoking. Having had BHR at a preceding visit was now borderline non-significant. A history of occupational sensitization had no statistically significant effect on ΔFEV1. However, when excluding duration of exposure, history of occupational sensitization predicted a statistically significant decrease of ΔFEV1 (data not shown). For all outcomes, we tested interactions between the different independent variables. None was statistically significant.

A c c e p t e d m a n u s c r i p t

DISCUSSION

In our study, BHR, respiratory symptoms and their simultaneous occurrence depended on duration of exposure, which confirms, for BHR, the cross -sectional results previously published on the inclusion data [START_REF] Bohadana | Respiratory symptoms and airway responsiveness in apparently healthy workers exposed to flour dust[END_REF]. Analysing respiratory symptoms , odds ratios for clinical atopy were high and statistically significant. FEV1

significantly decreased with duration of exposure and BHR at a preceding visit. This result persisted when accounting for the effect of BHR at the current visit. No effect of quantitative exposure level was observed.

In our study, the odds ratio of clinical atopy when analysing BHR was high although not significant. This is consistent with previous analyses of data of a cohort of apprentice bakers and pastry-makers, in which atopy or degree of sensitization was associated with incidence of BHR [START_REF] Tossa | Increase in exhaled nitric oxide is associated with bronchial hyperresponsiveness among apprentices[END_REF][START_REF] Wild | FeNO levels increase with degree of sensitisation in apprentices at risk of occupational asthma[END_REF]. Duration of follow-up in that study was shorter (18 months) than in the present study, and study populations were younger (mean age about 18 years) with no previous occupational exposure to flour. A possible explanation of the less clearcut effect of atopy could be that atopy has an effect on BHR in the first years of occupational exposure but no later effect on the onset of BHR. Indeed, the workers in our study had finished their apprenticeships since an average of 11 years at inclusion. However, it is noticeable that the prevalence of atopy was lower among bakers than in non-exposed subjects in our study. This suggests a selection effect with an early drop-out from the job of the subjects with atopy previous to the recruitment into the study. Such an effect was shown earlier in earlier studies among pastry-making or hairdressing apprentices [START_REF] Iwatsubo | Observatoire National des Asthmes Professionnels II. Bilan de la phase pilote et perspectives pour la surveillance des asthmes en lien avec le travail[END_REF][START_REF] Monso | Individual characteristics and quitting in apprentices exposed to high-molecular-weight agents[END_REF]. It could also be a healthy hire effect, by which subjects with childhood allergic diseases avoid high risk occupations, as suggested by results from a British cohort [START_REF] Butland | Job choice and the influence of prior asthma and hay fever[END_REF]. As in the Butland study, we used a definition of atopy based on self-reported childhood allergic diseases, which is different from the definition used in the former cited studies [START_REF] Monso | Individual characteristics and quitting in apprentices exposed to high-molecular-weight agents[END_REF][START_REF] Tossa | Increase in exhaled nitric oxide is associated with bronchial hyperresponsiveness among apprentices[END_REF][START_REF] Wild | FeNO levels increase with degree of sensitisation in apprentices at risk of occupational asthma[END_REF], based on specific IgE or SPT as in the standardized definition. This definition is the recommended one [START_REF] Johansson | A revised nomenclature for allergy. An EAACI position statement from the EAACI nomenclature task force[END_REF]), unfortunately we could not perform any of these investigations and used instead a definition based on a questionnaire. However, it has been suggested that a questionnaire may capture a more consistent chronic picture of atopy [START_REF] Hoppin | Questionnaire predictors of atopy in a US population sample: findings from the National Health and Nutrition Examination Survey, 2005-2006[END_REF] than these biochemical measures. Both prevalences of IgE sensitization or allergen skin test reactivity in (not elderly) adults have been found to decrease over 20 [START_REF] Amaral | Changes in IgE sensitization and total IgE levels over 20 years of follow-up[END_REF]10 (Warm et al. 2012) years of followup. However, self-reporting symptoms have other defaults such as recall bias.

A c c e p t e d m a n u s c r i p t

In our study, we failed to observe a statistically significant effect of occupational sensitization on BHR, respiratory symptoms or ∆FEV1, although in the literature, sensitization to occupational allergens has been observed to be a risk factor for work-related symptoms in bakers [START_REF] Droste | Occupational exposure among bakery workers: impact on the occurrence of work -related symptoms as compared with allergic characteristics[END_REF][START_REF] Harris-Roberts | Sensitization to wheat flour and enzymes and associated respiratory symptoms in British bakers[END_REF][START_REF] Hur | Prevalence of work-related symptoms and serum-specific antibodies to wheat flour in exposed workers in the bakery industry[END_REF][START_REF] Jacobs | Wheat allergen exposure and the prevalence of work-related sensitization and allergy in bakery workers[END_REF][START_REF] Ngahane | Prevalence and risk factors for allergic rhinitis in bakers in Douala[END_REF][START_REF] Olivieri | Wheat IgE profiling and wheat IgE levels in bakers with allergic occupational phenotypes[END_REF] or for exercise-induced bronchoconstriction [START_REF] Minov | Exercise -related respiratory symptoms and exercise-induced bronchoconstriction in industrial bakers[END_REF]. Moreover, high levels of flour-specific IgE have been found to be good predictors of a positive specific inhalation challenge test [START_REF] Van Kampen | Prediction of challenge test results by flour-specific IgE and skin prick test in symptomatic bakers[END_REF]. A cross-sectional study in bakers did not observe a relationship between occupational sensitization and BHR [START_REF] Storaas | Bronchial responsiveness in bakery workers: relation to airway symptoms, IgE sensitization, nasal indices of inflammation, flour dust exposure and smoking[END_REF]). The authors hypothesized a non-IgE mediated mechanism causing BHR and airway symptoms, based on their former works [START_REF] Storaas | Occupational rhinitis: diagnostic criteria, relation to lower airway symptoms and IgE sensitization in bakery workers[END_REF].

Concerning our study, we could only test three occupational allergens (wheat flour, rye flour and alpha amylase)

and might have missed other relevant ones. As mentioned in the introduction, the use of a large variety of additives has developed during the last decades [START_REF] Quirce | Diagnosis and management of grain-induced asthma[END_REF]. However, this does not mean there was no effect of occupational sensitization in our study. Indeed, when no exposure variables were present in the model (data not shown), FEV1 was significantly lower among workers with occupational sensitization. The fact that occupational sensitization was only a predictor of a decline in FEV1 when duration of exposure was not included might mean that exposure predicts this decline and sensitization independently but that the latter is not an intermediate in the causal relation between exposure and decline in FEV1. In our study, the prevalence of BHR and symptoms depended on gender and atopic background. Similar to other studies, a higher prevalence of BHR has been observed in women than in men [START_REF] Brutsche | Bronchial hyperresponsiveness and the development of asthma and COPD in asymptomatic individuals: SAPALDIA cohort study[END_REF] and in atopics than in non-atopics in a previous study in bakers [START_REF] Pavlovic | Bronchial hyperactivity in bakers and its relation to atopy and skin reactivity[END_REF]).

We identified only few studies exploring the longitudinal decline of FEV1 according to BHR [START_REF] Abramson | Predictive value of non-specific bronchial challenge testing for respiratory symptoms and lung function in aluminium smelter workers[END_REF][START_REF] Aldrich | Bronchial Reactivity and Lung Function After[END_REF][START_REF] Brutsche | Bronchial hyperresponsiveness and the development of asthma and COPD in asymptomatic individuals: SAPALDIA cohort study[END_REF]. A population-based study [START_REF] Brutsche | Bronchial hyperresponsiveness and the development of asthma and COPD in asymptomatic individuals: SAPALDIA cohort study[END_REF]) found a relationship between FEV1 and BHR but did not separate past and concurrent BHR as determinant s of FEV1. In new-employed aluminium smelter workers, Abramson et al. did not find any predictive value of BHR at start of employment concerning the decline of lung function [START_REF] Abramson | Predictive value of non-specific bronchial challenge testing for respiratory symptoms and lung function in aluminium smelter workers[END_REF]. The median of follow-up was 4 years, which might be too short to be able to observe a detectable functional decline. A study with a longer follow-up in firefighters [START_REF] Aldrich | Bronchial Reactivity and Lung Function After[END_REF] found that BHR at 12-year-follow-up was associated with an estimated 15.4 ml/y greater FEV1 decline than experienced in case of no BHR at follow-up. In our study, the sum of the pre-visit and current BHR effects corresponds to a 180 ml approximated FEV1 loss during 13 years of follow-up, equivalent to a 14 ml/year loss. This is roughly in line with Aldrich et al.'s results [START_REF] Aldrich | Bronchial Reactivity and Lung Function After[END_REF], although the context of exposure is very different: a massive exposure to inorganic dust and gases for the World Trade Center's collapseexposed firefighters and a chronic exposure to mostly organic dust for the bakers. However, the analyses in firefighters were adjusted on the changes in FEV1 shortly before and after the 9/11/2001, controlling for the effect of this massive exposure. In our study, duration of exposure predicted an additional decrease of near 20 ml/y in FEV1, independently of the effect of BHR. In Aldrich et al.'s study, the exposure after the World Trade Center's collapse was not taken into account, 76% of the participants being retired at the end of follow-up [START_REF] Aldrich | Bronchial Reactivity and Lung Function After[END_REF]).

In our study, we did not observe any dose-response relationships using the quantitative exposure. As the exposure estimation was based on on-going exposure measurements and assessments by industrial hygienists, it is unlikely that it was due to exposure misclassification. A more likely mechanism would be that the high -exposed tasks are only performed by workers, with potentially low sensitization, who withstand these high exposure levels. It could be that workers themselves move to lower-exposed tasks or are more careful, e.g. wearing protective equipment from the appearance of the first work-related symptoms.

Our study has several strengths. First, we followed up the population over 13 years with the same study staff using the same protocol. Second, for ten bakers (9%) only, the reason for drop-out was unknown, and two bakers only dropped out for reasons related to occupational health . This loss to follow-up rate is quite low when compared to studies on apprentices (e.g. 20% in [START_REF] Tossa | Increase in exhaled nitric oxide is associated with bronchial hyperresponsiveness among apprentices[END_REF]). This could be due to the fact that all the data collection was done at the workplace. The corresponding drawback was that we had to use an abbreviated version of the MBC test for ethical reason. Subjects lost to follow-up in our study were older, had a shorter duration of employment but had a greater FEV1 than the other bakers (data not shown). However, no statistical difference between them and the other bakers was observed. Third, an attrition bias during the follow-up seems unlikely as no difference in health respiratory status was found at the last visit.

The main weakness of this study is the relatively low sample size, due to small-medium size of the working resources in that activity sector in France. In particular the low number of subjects in the non-exposed might be considered a problem. One must however bear in mind that the results of the present study do not rely on an exposed vs. non-exposed comparison. The non-exposed population is taken as a baseline exposure category in our search of dose-response analyses, in the present paper with duration of exposure and levels of exposure . Thus the inclusion of the non-exposed subjects strengthened the power of detecting such a dose-response. Moreover, the A c c e p t e d m a n u s c r i p t inclusion allowed to better take into account confounders. Indeed (data not shown), repeating the analyses without the non-exposed confirmed the effect of the duration of exposure. On the other hand, some confounding factors (pack-years, history of BHR) were no longer significant.

A second possible weakness is the absence of systematic individual exposure measurements which led us to use a semi-quantitative exposure levels in three categories in whose reliability we could trust. The sketchy individual exposure measurements would have been unreliable. Moreover the measurements were obtained using a 1L/min flowrate, which were up to date at the time of the measurements but are lower than the flowrates (2.5-3L/ min ) which are recommended nowadays. The precision of these suboptimal measurements was however deemed sufficient for the classification in the three exposure categories used.

Finally, a selection effect with an avoidance of the job or an early drop-out from the job of the subjects with atopy cannot be excluded, as discussed above. Despite this healthy worker selection effect, we could observe detrimental effects of exposure on bronchial responsiveness and on FEV1.

CONCLUSION

In flour-exposed industrial bakers, length of exposure and smoking are long-term determinants of BHR and of the decrease in FEV1. BHR at a preceding visit predicted lower FEV1 even when accounting for the effect of BHR at the current visit. Prevention is still needed for these workers.

ETHICS APPROVAL: The study was approved by the relevant French ethical committee (Comité consultatif de protection des personnes dans la recherche biomédicale de Lorraine) and the French data protection agency (Commission nationale de l'informatique et des libertés) and has therefore been performed in accordance with the ethical standards laid down in the 1964 Declaration of Helsinki and its later amendments . All the subjects gave free written and informed consent to participate.

Figure 1

 1 Figure 1 Flow chart of the study. *drop-out unrelated to health A c c e p t e d m a n u s c r i p t

  A c c e p t e d m a n u s c r i p t were adjusted for the technician performing the spirometry. Age, gender and height were taken into account by subtracting the predicted FEV1 from the observed values.

  A c c e p t e d m a n u s c r i p t

Figure caption Fig 1 :

 caption1 Figure caption

Table 1 .

 1 Basic description of population study

	Exposed subjects	Non-exposed

Table 2 .

 2 Description of the health outcomes

		Bronchial	Respiratory	FEV1 in % predicted c	Total
		hyperresponsiveness a	symptoms b	(mean [SD])
	At inclusion
	Non-exposed workers	3 (7.9%)	13 (34.2%)	104.5 [15.0]	38
	Missing data	2 (5.3%)	0	0
	Exposed workers	26 (23.6%)	42 (38.2%)	99.6 [13.8]	110
	Exposure level (mg/m 3 )
	Missing data	4 (3.6%)	0	0
	<1	5 (18.5%)	5 (18.5%)	95.6 [15.9]	27
	[1-5]	19 (25.7%)	36 (48.7%)	101.1 [13.3]	74
	>5	2 (22.2%)	1 (11.1%)	99.7 [9.9]	9
	Duration of exposure (years)
	<10	10 (18.9%)	19 (35.9%)	101.5 [12.3]	53
	[10-20]	8 (19.5%)	14 (34.1%)	100.2 [15.1]	41
	>20	8 (50.0%)	9 (56.2%)	91.9 [13.3]	16
	At the last visit
	Non-exposed workers	4 (10.5%)	15 (39.5%)	101.5 [14.9]	38
	Missing data	5 (13.2%)	0	1 (2.6%)
	Exposed workers	40 (36.4%)	57 (51.8%)	99.8 [15.7]	110
	Exposure level (mg/m 3 )
	Missing data	0	0	0
	<1	11 (44.0%)	12 (48.0%)	101.9 [19.0]	25
	[1-5]	27 (34.2%)	44 (55.7%)	99.2 [15.0]	79
	>5	2 (33.3%)	1 (16.7%)	98.6 [11.7]	6
	Duration of exposure (years)
	<10	10 (31.2%)	16 (50.0%)	104.7 [13.0]	32
	[10-20]	13 (32.5%)	21 (52.5%)	101.5 [15.4]	40
	>20	17 (44.7%)	20 (52.6%)	93.9 [16.6]	38
	a

Table 3 .

 3 Models for the health outcomes according to a mixed logistic regression (OR [95%CI])

		BHR	Respiratory	BHR and respiratory
			symptoms a	symptoms
		n = 103	n = 198 respiratory	n = 51 BHR and respiratory
		BHR/445 visit-	symptoms/462 visit-	symptoms/445 visit-subjects
		subjects	subjects
	Gender: women vs men	13.1 [1.9;91.8]	2.2 [0.3;18.4]	2.9 [0.5;17.1]
	Cumulative smoking (10 pack-	1.5 [1.0;2.4]	1.2 [0.7;1.9]	1.3 [0.9;1.9]
	years)
	Clinical atopy: yes vs no	5.2 [0.8;31.9]	16.5 [1.7;157.1]	5.8 [1.1;31.4]
	History of occupational	1.4 [0.4;4.7]	1.3 [0.3;5.1]	1.7 [0.5;6.2]
	sensitization b : yes vs no
	Duration of exposure (decade)	2.2 [1.0;5.0]	3.6 [1.4;9.4]	3.8 [1.5;9.5]
	Mean exposure level (mg/m 3 )
	Reference (non-exposed)	1	1	1
	<1	7.9 [0.7;92.3]	0.6 [0.1;5.3]	0.7 [0.1;8.2]
	[1-5]	4.5 [0.4;47.2]	0.9 [0.1;6.2]	0.5 [0.0;4.9]
	>5	2.6 [0.2;33.8]	0.2 [0.0;2.2]	0.5 [0.0;6.3]
	BHR : bronchial hyperresponsiveness defined as fall in FEV1 of 20% or more during the MBC test

a : respiratory symptoms defined as a wheezing chest or a runny nose in absence of any cold b : sensitization to occupational allergens at least once during follow-up

Table 4 .

 4 Models for ΔFEV1 (Observed-predicted according to gender, height and age) in ml according to a multiple linear mixed analysis (β coefficient [95%CI]).

		Model 1	Model 2
		N= 466 visits for 146 subjects	N= 445 visits for 146
			subjects
	Cumulative smoking (10 pack-years)	-61.8 [-111.8 , -11.7]	-57.9 [-108.3 , -7.4]
	Pre-visit BHR: yes vs no	-81.4 [-161.8 , -0.1]	-76.5 [-154.8 , 1.8]
	Having BHR at the current visit: yes vs no	-95.3[-168.0 , -22.7]	-93.6 [-163.9 , -23.2]
	History of occupational sensitization a : yes vs no		-63.2 [-165.8 , 39.4]
	Duration of exposure (decade)	-193.5 [-268.0 , -119.1]	-188.2 [-254.1 , -122.2]
	Mean exposure level (mg/m 3 )
	Reference (non-exposed)	1
	<1	96.8 [-122.4 , 316.0]
	[1-5]	94.1 [-120.4 , 308.6]
	>5	14.8 [-74.0 , 370.7]
	CI: Confidence Interval
	β coefficient: multiple linear regression coefficient
	BHR: bronchial hyperresponsiveness defined as fall in FEV1 of 20% or more during the MBC test

a : sensitization to occupational allergens at least once during follow-up

Acknowledgments

The authors are grateful to all workers who participated in th e study. They thank all company directors, managers and occupational physicians who were involved. They also thank Régis Colin from INRS for his precious help in the data management.

Declarations


D Vola

M Raous

J A C Martins

Friction and instability of steady sliding: squeal of a rubber/glass contact

The dynamic instability of steady sliding states of finite dimensional frictional contact systems with non-linear elastic behavior is analyzed. An algorithm for the computation of those steady sliding states and a sufficient condition for their instability, based on the resolution of a generalized eigenvalue problem, are presented. Flutter instabilities due to the non-associative character of the Coulomb friction law are shown to occur for a finite element model of a rubberlike waist seal sliding on a glass window that is known to generate squeal noise. The consequences of those flutter instabilities are assessed by computing various finite element dynamic solutions in the neighborhood of steady sliding.

INTRODUCTION

This work is part of a research effort devoted to the study of instability phenomena in frictional contact problems. More specifically, it deals with the generation of noise during the sliding of a deformable body on a rigid obstacle. The case of a rubber/glass contact between a waist seal and a car window is discussed.

In a previous work [START_REF] Martins | Dynamic stability of finite dimensional linearly elastic systems with unilateral contact and Coulomb friction[END_REF] , the dynamic stability of a finite dimensional frictional contact system with linear elastic behavior was discussed. The theory was applied to study the possible origines for the occurrence of jumps on the time evolution of the friction force and the propagation of stress waves along the contact zone, in the course of the sliding of a polyurethane block on a flat araldite plate [START_REF] Martins | Dynamic stability of finite dimensional linearly elastic systems with unilateral contact and Coulomb friction[END_REF][START_REF] Raous | Stress waves in a sliding contact. Part 1 : modelling[END_REF][START_REF] Vola | Bifurcations and instabilities in some finite dimensional frictional contact problems[END_REF][START_REF] Zeghloul | Stress waves in a sliding contact. Part 1 : experimental study[END_REF] , . A generalization of this work to finite dimensional nonlinear elastic systems in frictional contact with curved obstacles may be found in [START_REF] Martins | Stability of finite dimensional systems with unilateral contact and friction: non-linear elastic behaviour and obstacle curvature[END_REF] .

In the present work, a nonlinear material behavior is considered (finite elasticity, Mooney Rivlin strain energy), and stability analyses are performed for steady sliding equilibrium states. A fixed point algorithm on the tangential contact stresses is used to compute the steady sliding solutions in a range of values of the coefficient of friction, which is the essential parameter in the stability analyses. The occurrence of dynamically unstable steady sliding states (existence of smooth growing dynamic solutions with perturbed initial conditions arbitrarily close to steady sliding) is predicted by solving appropriate generalized eigenvalue problems. The two matrices in these problems are finite element matrices depending on the frictional contact conditions. When a flutter instability is detected, its consequences are studied by performing dynamic analyses of the system behavior. The formulation of the finite dimensional dynamic problem is done in terms of differential measures, leading to a numerical method that handles the velocity jumps on the solution.

Finite dimensional stability analyses of steady sliding states were also performed in [START_REF] Martins | A study of static and kinetic friction[END_REF][START_REF] Oden | Models and computational methods for dynamic friction phenomena[END_REF] , for some metallic bodies, under the assumption of rigid or linear elastic behavior, and adopting a nonlinear compliance law for the normal behavior of the interface. Although many aspects of the present work are quite distinct from those of [START_REF] Martins | A study of static and kinetic friction[END_REF][START_REF] Oden | Models and computational methods for dynamic friction phenomena[END_REF] (geometry and boundary conditions of the bodies, material and interface properties, and also the computational algorithms) it is clear that the cause for the flutter instabilities detected in [START_REF] Martins | A study of static and kinetic friction[END_REF][START_REF] Oden | Models and computational methods for dynamic friction phenomena[END_REF] and in the present paper is the same : the non-symmetry (non-self-adjointness) of Coulomb's friction law. Furthermore, despite all those distinctions, the consequences of the flutter instabilities found in the present paper have some qualitative features that clearly resemble some of those discussed in [START_REF] Martins | A study of static and kinetic friction[END_REF] .

In section 2, the continuum formulation of the dynamic problem in finite elasticity with frictional unilateral contact is presented. Section 3 is devoted to formulate and approximate the finite dimensional dynamical problem. Section 4 focuses on the direct computation of the steady sliding solution, and the stability analysis is presented in section 5. In section 6, an application to the waist seal/car window squeal (a rubber/glass contact) is presented. Flutter instabilities with frequencies in the acoustic range, may indeed occur in the studied system.

THE FRICTIONAL CONTACT PROBLEM IN FINITE ELASTICITY

Finite elasticity for rubberlike materials

Having in mind the application to a rubberlike waist seal and neglecting the viscosity, an incompressible hyperelastic isotropic material is considered. Let B be a deformable body occupying the domain Ω at a reference configuration C. The motion of the body B is given by the map ϕ

ϕ ϕ ϕ : Ω × IR → IR 3 , (X, τ ) → x = ϕ ϕ ϕ ϕ(X, τ ) = ϕ ϕ ϕ ϕ τ (X)
. The material constraint of incompressibility is classically written J = det(F) = 1, where F = ∇ X (ϕ ϕ ϕ ϕ τ ) is the deformation gradient and det denotes the determinant of a second order tensor. In order to overcome locking effects, we choose, among the numerous formulations in the literature [START_REF] Gadala | Numerical solutions of nonlinear problems of continua -II. A survey of incompressibility constraints and software aspects[END_REF] , the three-field formulation for quasi-incompressibility proposed by J.C.Simo and R.L.Taylor 24 . This formulation was also used by G.Björkman et al. [START_REF] Björkman | Sequential quadratic programming for non-linear elastic contact problems[END_REF] for quasistatic frictionless contact problems and is now applied to the present dynamical problem with frictional contact. It is based on the split proposed by P.J.Flory 4 for the deformation gradient F and, consequently, for the (right) Cauchy-Green tensor C, into isochoric and volumetric parts :

F = J 1/3 F C = J 2/3 C Ĩ3 = 1, (1)
where Ĩ3 is the third principal invariant of C. The volume change is characterized by using a separate field θ, and a third unknown field (a Lagrange multiplier p) enforces the constraint :

θ = J. (2)
Moreover the volumetric effects on the strain energy and the deviatoric ones are supposed to be uncoupled. This leads to the following expression for the strain energy :

W ( Ĩ1 , Ĩ2 , θ) = W ( Ĩ1 , Ĩ2 ) + U (θ), (3)
where Ĩ1 and Ĩ2 are the first and second principal invariants of C. We have chosen the Mooney Rivlin strain energy that is sufficient for studying the waist seal deformation (the maximum strain numerically obtained for the waist seal is of 50%) :

W = a 1 ( Ĩ1 -3) + a 2 ( Ĩ2 -3), (4)
where a 1 and a 2 are the material coefficients. Then the constitutive behavior of the material is characterized by :

S = pJC -1 + 2J -2/3 Dev( ∂ W ∂ C ), (5)
p = U ′ (θ), (6)
and θ satisfies the constraint (2). S is the second Piola-Kirshoff stress tensor,

Dev(A) = A -1 3 (A : C)C -1 and ∂ W /∂ C = (a 1 + a 2 Ĩ1 )ll -a 2 C. Γ u t V obs τ ϕ n ϕ ( τ τ Γ Reference configuration Current configuration ϕ ( Ω) Γ ) c Γ c Γ s N Ω g Figure 1:
Deformable body in contact with a flat obstacle.

The volumetric part U is similar to a penalty function. It should verify the following properties : U is non-negative over IR+, U (1) = 0 and lim J→0 U (θ) = ∞. That is the case with the following choice (κ is the penalty parameter homogeneous to a bulk modulus) :

U (θ) = κ[ 1 2 (θ 2 -1) -ln(θ)]. (7)

Frictional contact conditions

We restrict this study to the two-dimensional problem of the contact between a deformable body and a flat and rigid obstacle that moves with a constant velocity v obs = v obs t along a direction t of its surface (see figure 1). We denote by Γ c the part of the boundary of the reference configuration of the body that contains the particles that may eventually establish contact with the obstacle. The gap between the obstacle surface and the current position x of such contact candidate particles (x = ϕ ϕ ϕ ϕ(X, τ ) for X belonging to Γ c ) is represented by G(X, τ ) or, equivalently, by g(x, τ ) = G(ϕ ϕ ϕ ϕ -1 (x, τ ), τ ). The component along the direction tangent to the obstacle of the velocity of a contact candidate particle relative to the obstacle is represented by

V (X, τ ) = φ ϕ ϕ ϕ(X, τ ).t -v obs or, equivalently, by v(x, τ ) = V (ϕ ϕ ϕ ϕ -1 (x, τ ), τ ). The contact reactions r on ϕ ϕ ϕ ϕ τ (Γ c ) are related to the Cauchy stress tensor σ σ σ σ = 1/JFSF T by r = σ σ σ σn Γ , (8)
where n Γ is the normal to the boundary of the body in its current configuration, and the decomposition in normal and tangential components to the obstacle is performed : r = r n n + r t t, r n = r.n and r t = r.t. The contact is supposed to be unilateral :

g ≥ 0, r n ≥ 0 and r n g = 0, (9)
and the Coulomb law of friction is adopted :

|r t | -µr n ≤ 0, µr n |v| = -r t v, ( 10
)
where µ is the friction coefficient. Unlike most of the classical analyses that associate friction induced instabilities to variations of the coefficient of friction, µ is taken as a r (R) Cauchy (nominal) vector of reaction stresses on the contact surface, g (G)

Eulerian (Lagrangean) descriptions of the current gap, v obs Absolute value of the tangential velocity of the flat obstacle, v (V ) Eulerian (Lagrangean) descriptions of the current sliding velocity, P Vector of dimension 2N c containing the nodal contact reaction forces for the finite element model, G

Vector of dimension N c containing the gaps between the contact candidate nodes and the flat obstacle for the reference configuration of the f.e. mesh, H n (H t ) Transformation matrices of dimension N c × N that relate the normal (tangential) displacements of the contact candidate nodes of the f.e. mesh to the global vector of nodal displacements.

Table 1: Nomenclature for the main contact variables constant in the present work. Using the framework of convex analysis 7 , the friction condition (10) becomes :

-r t ∈ ∂φ rn (v), (11)
where φ rn (y) = µr n |y| and ∂φ rn denotes the sub-differential of φ rn . The Coulomb friction condition may also be written as a maximal dissipation principle :

|r t | ≤ µr n and ∀r * t such that |r * t | ≤ µr n , (r * t -r t )v ≥ 0. ( 12
)
Both contact and friction conditions can be expressed in the reference configuration by using the nominal contact reactions :

R = FS N Γ , (13)
where N Γ is the normal to the boundary Γ of the reference configuration of the body, and the decomposition R = R n n + R t t, R n = R.n and R t = R.t can also be used :

G ≥ 0, R n ≥ 0 and R n G = 0, ( 14
)
|R t | ≤ µR n and ∀R * t such that |R * t | ≤ µR n , (R * t -R t )V ≥ 0. ( 15
)
A nomenclature for the main contact variables can be found in table 1.

THE DYNAMIC SOLUTION

Formulation and semi-discretization

In this subsection, we are interested in the computation of the dynamic behavior of the solid B. The perfect unilateral contact and Coulomb friction conditions induce mathematical difficulties that are still unresolved for continuum bodies. Thus the equations of elastodynamics and the corresponding variational problem are formally stated for sufficiently regular solutions. After space discretization we will state a problem that takes into account discontinuities of the contact velocity [START_REF] Jean | Dynamics in the presence of unilateral contact and dry friction: a numerical approach[END_REF] . At any instant τ , the equation of linear momentum balance is :

ρ φ ϕ ϕ ϕ -Div(FS) = f v , ( 16
)
where ρ is the density of the body in its reference configuration and f v is the external force that acts on Ω. The boundary conditions are given by the prescribed configuration on

Γ u ϕ ϕ ϕ ϕ τ = ϕ ϕ ϕ ϕ imp on Γ u , (17)
by the surface tractions f s on Γ s ,

(FS)N Γ = f s , (18)
where N Γ is the outward unit normal vector to Γ, and by the frictional contact conditions on Γ c (Γ = Γu ∪ Γs ∪ Γc and Γ u , Γ s and Γ c being disjoint). The initial conditions are given by the initial positions and the initial velocities.

Weak formulation

Considering the specific form (7) of the volumetric part U of the strain energy W , we suppose that no violation of the kinematic constraint det(∇ X (ϕ ϕ ϕ ϕ τ )) > 0 is possible [START_REF] He | Unilateral contact of two solids subject to large deformations: formulation and existence results[END_REF] . The set of admissible configurations is simply :

W = {ϕ ϕ ϕ ϕ τ : Ω -→ IR 2 such that ϕ ϕ ϕ ϕ τ (X) = ϕ ϕ ϕ ϕ imp ∀X ∈ Γ u }. ( 19
)
The set of admissible variations to the reference configuration is :

T 0 W = { η η η η : Ω -→ IR 2 such that η η η η(X) = 0 ∀X ∈ Γ u }. (20)
The set of admissible hydrostatic pressures and volume changes is :

V = {θ * : Ω -→ IR}. (21)
Finally, the set of admissible normal contact stresses is :

C n = {R * n : Γ c -→ IR such that R * n ≥ 0}. ( 22
)
Except for the inertia forces, all the non-dissipative effects can be put together in the functional Π proposed by G.Björkmann et al. [START_REF] Björkman | Sequential quadratic programming for non-linear elastic contact problems[END_REF] :

Π(ϕ ϕ ϕ ϕ τ , θ, p, R n ) = Π int (ϕ ϕ ϕ ϕ τ , θ, p) -Π ext (ϕ ϕ ϕ ϕ τ ) - Γc R n G dΓ, (23) with
(ϕ ϕ ϕ ϕ τ , θ, p, R n ) ∈ W × V × V × C n ; Π int denotes an augmented potential of the internal forces : Π int (ϕ ϕ ϕ ϕ τ , θ, p) = Ω W ( Ĩ1 , Ĩ2 ) + U (θ) + p(J -θ) dV , (24)
and Π ext denotes the potential of the external forces :

Π ext (ϕ ϕ ϕ ϕ τ ) = Ω f v ϕ ϕ ϕ ϕ τ dV + Γs f s ϕ ϕ ϕ ϕ τ dΓ. (25)
The following pseudo-potential represents the dissipation due to friction :

D(R n , φ ϕ ϕ ϕ τ -Φ Φ) = Γc µR n |( φ ϕ ϕ ϕ τ -Φ Φ).t| dΓ, (26)
where Φ Φ is an application that verifies Φ Φ ∈ W and Φ Φ(X) = v obs on Γ c . Taking into account the Coulomb friction condition (11) and using the notion of sub-differential of a convex function (see [START_REF] Rockafellar | Convex analysis[END_REF] ), the variational problem can be written as an inclusion 7 :

-Π ,ϕ (ϕ ϕ ϕ ϕ τ , θ, p, R n ) -F ine ∈ ∂ 2 D(R n , φ ϕ ϕ ϕ τ -Φ Φ), (27)
together with the following variational equations and inequality :

∀θ * ∈ V, Π int,θ θ * = Ω (U´(θ) -p)θ * dV = 0, (28)
∀p * ∈ V, Π int,p p * = Ω (J -θ)p * dV = 0, (29)
∀R * n ∈ C n , Γc (R * n -R n )G dΓ ≥ 0, ( 30
)
where F ine denotes the inertia forces and ∂ 2 D the sub-differential of D with respect to its second argument. Using the definition of the sub-differential of a convex function, the inclusion ( 27) may be replaced by the variational inequality :

∀ϕ ϕ ϕ ϕ * ∈ W, -< Π ,ϕ + F ine ; ϕ ϕ ϕ ϕ * -φ ϕ ϕ ϕ τ >≤ D(R n , ϕ ϕ ϕ ϕ * -Φ Φ) -D(R n , φ ϕ ϕ ϕ τ -Φ Φ). (31)
If the Coulomb friction condition is expressed in terms of the maximal dissipation principle (12) and tangential nominal contact reactions R t are considered as an additional unknown field, then the variational inequality (31) is replaced by the following variational equation (d'Alembert's principle) :

∀ η η η η ∈ W, (Π int,ϕ -Π ext,ϕ + F ine ). η η η η - Γc (R n n + R t t). η η η η dΓ = 0, (32)
together with the variational inequality (Coulomb friction condition) :

∀R * t ∈ C t (R n ), Γc (R * t -R t )V dΓ ≥ 0, ( 33
) with Π int,ϕ . η η η η = Ω FS : ( ∂ η η η η ∂X )dV , (34)
Π ext,ϕ . η η η η = Ω f v η η η η dV + Γs f s η η η η dΓ, (35)
and where C t (R n ) is the set of admissible tangential contact reactions :

C t (R n ) = {R * t : Γ c -→ IR such that |R * t | ≤ µR n }. (36)
Summarizing, the dynamical frictional contact problem P1 is governed at each time τ by the set of variational equations and inequalities ( 28)-( 31) involving the unknown fields ϕ ϕ ϕ ϕ, θ, p and R n or, equivalently, by the variational statements ( 28)-( 30), ( 32) and (33) involving the unknown fields ϕ ϕ ϕ ϕ, θ, p, R n and R t .

Space discretization

Many details about the space dicretization associated to the three-field formulation may be found in [START_REF] Björkman | Sequential quadratic programming for non-linear elastic contact problems[END_REF][START_REF] Liu | 3d finite element analysis of rubber-like materials at finite strains[END_REF][START_REF] Simo | Quasi-incompressible finite elasticity in principal stretches. Continuum basis and numerical algorithms[END_REF] . The number of unknowns may be reduced by eliminating the two auxilary fields θ and p at the element level. Indeed the most important feature is that the shape functions of θ and p are chosen to be discontinuous over the element boundaries. That is the case of the finite element Q4P0 used in this study (four nodes for the displacement interpolation and one node for the interpolation of the volume change and the hydrostatic pressure). If θ h e and p h e denote the discretized forms of the volume change and of the pressure on some finite element Ω h e , the equations ( 28) and (29) yield the following relations :

θ h e = V ol(ϕ ϕ ϕ ϕ τ (Ω h e )) V ol(Ω h e ) , p h e = U´(θ h e ). ( 37
)
Then, at the element level, θ h e and p h e are replaced in the equations of motion by the relations (37) [START_REF] Simo | Quasi-incompressible finite elasticity in principal stretches. Continuum basis and numerical algorithms[END_REF] .

The integrals over the boundary Γ are discretized by using the trapezoidal rule. The integration points are the nodes of the displacement interpolation. We denote by P n the vector (of dimension N c ) of normal nodal contact forces (respectively P t for the tangential nodal contact forces) where N c is the number of contact candidate nodes. The vector P of dimension 2 * N c puts together the components of P n and P t . We denote by U the vector of dimension N of nodal displacements (N is the total number of degrees of freedom). The normal and tangential components of the displacement of the contact node number i are given by :

U n i = H n i U U t i = H t i U, (38)
where

H
C h n = {P n ∈ IR Nc such that P n i ≥ 0 ∀i = 1, .., N c }, C h t (P n ) = {P t ∈ IR Nc such that |P t i | ≤ µP n i ∀i = 1, .., N c }. (39)
Using the previous notations, the discretized form of the variational inequality (30) of the unilateral contact condition is :

∀P * n ∈ C h n (P * n -P n ) T (H n U + G) ≥ 0, ( 40
)
where each component of the vector G (of dimension N c ) is equal to the corresponding reference gap. The discretized form of the variational inequality (33) of the friction conditions is :

∀P * t ∈ C h t (P n ) (P * t -P t ) T (H t U -V obs ) ≥ 0, ( 41
)
where the components of the vector V obs (of dimension N c ) are equal to v obs . The semidiscrete equation of motion is :

M Ü + F int (U) -F ext = H T P (42)
where M, F int and F ext denote the consistent mass matrix, the internal forces and the external forces, respectively. The semi-discrete Problem P2 is thus given by the equation of motion (42), the initial displacements U 0 and the initial velocities U0 together with the frictional contact conditions (40) and (41).

Finite dimensional dynamical contact problems are not regular because of the discontinuities of the gap velocity when impact occurs. In [START_REF] Jean | Dynamics in the presence of unilateral contact and dry friction: a numerical approach[END_REF][START_REF] Jean | The non smooth contact dynamics method[END_REF] the problem is set so as to take into account those discontinuities. Following J.J.Moreau 7 , we make the assumption that the velocity has a variation bounded in time so as to define at any instant a left Uand a right U+ velocity. We suppose also that the velocity is right-continuous and we denote by U = U+ the velocity at any instant. The acceleration should be understood as a differential measure d U :

]a,b] d U = U+ (b) - U+ (a) (43)
and the equation of motion (42) as a measure differential equation :

Md U + F int (U)dτ -F ext dτ -H T Pdν = 0 ( 44
)
where dν is some non-negative measure relative to which d U possesses a density function and P is the representative of the local density of impulses exerted when contact occurs.

Approximation in time

Following S.Pandit and S.Deo 18 , one obtains :

       M( U(τ ) -U(0)) = τ 0 (F ext -F int (U))dτ + ]0,τ ] H T Pdν U(τ ) -U(0) = τ 0 Udτ (45)
After dividing [0, T ] into sub-intervals, several time approximations of (44) may be used, see [START_REF] Vola | Consistent time discretization for a dynamical frictional contact problem and complementarity techniques[END_REF] . We have chosen the following one :

                   τ m+1 τ m (F ext -F int (U ))dτ ≈ ∆t [θ 1 (F m+1 ext -F int (U m+1 )) +(1 -θ 1 )(F m ext -F int (U m ))] τ m+1 τ m Udτ ≈ ∆t[θ 2 Um+1 + (1 -θ 2 ) Um ] ]τ m ,τ m+1 ] H T Pdν ≈ ∆tH T Pm+1 ( 46
)
where 2 m denotes the approximation of 2(τ m ), ∆t = τ m+1 -τ m , θ 1 and θ 2 are two numerical coefficients in the interval [0,1] and Pm+1 is the mean value of the impulse density in the time interval.

Using the second approximation of (46), the contact condition is written with the velocity as unknown :

     Pm+1 n ∈ C h n and ∀P * n ∈ C h n (P * n - Pm+1 n ) T (H n Um+1 + Gm ) ≥ 0 Gm = (1 -1 θ 2 )H n Um + 1 ∆tθ 2 (H n U m + G) (47)
The discretized Coulomb condition involves the tangential relative velocity H t Um+1 -V obs and the average contact force Pm+1 . Finally the discrete non-linear problem is the following one :

Problem P3 : find Um+1 ∈ IR N , Pm+1 n ∈ C h n and Pm+1 t ∈ C h t ( Pm+1
n
) such that :

     M Um+1 + ∆tθ 1 (F int (l( Um+1 )) -F m+1 ext ) = ∆tF m res + ∆tH T Pm+1 ∀P * n ∈ C h n (P * n - Pm+1 n ) T (H n Um+1 + Gm ) ≥ 0 ∀P * t ∈ C h t ( Pm+1 n ) (P * t - Pm+1 t ) T (H t Um+1 -V obs ) ≥ 0 (48) with l( Um+1 ) = 1 ∆tθ 2 (U m+1 -U m ) -1-θ 2 θ 2 Um F m res = 1 ∆t M Um + (1 -θ 1 )[F m ext -F int (l( Um ))]. ( 49
)

Resolution

In order to linearize the problem (48) which is constrained by variational inequalities, it is formulated in terms of generalized equations. The variational inequality (C is a convex set) :

z ∈ C and ∀z * ∈ C F(z)(z * -z) ≥ 0 ( 50
)
is equivalent to the generalized equation :

-F(z) ∈ ∂ψ C (z) (51)
where ψ C is the indicator function of the convex set C. Following A.Klarbring and G.Björkmann 11 , we suppose that implicit sets such as

C h t ( Pm+1
n
) may also be considered. For simplicity, the superscript m+1 is now omitted. In the case of the dynamical problem, we consider :

C = IR N × C h n × C h t ( Pn ) and z =   U Pn Pt   , F(z) =   L( U, Pn , Pt ) Gm n + H n U V obs -H t U   , (52)
with

L( U, Pn , Pt ) = F m res + H T n Pn + H T t Pt -1 ∆t M Um -θ 1 F int (l( U)).
The Newton-Raphson algorithm may be applied to the generalized equations [START_REF] Cottle | The linear complementarity problem[END_REF] .

Considering the properties of H, the linearisation leads to the following problem.

Problem P4 : Assuming that the evolution is known until the instant τ m and that ( U0 , ( Pn ) 0 , ( Pt ) 0 ), .., ( Uk , ( Pn

) k , ( Pt ) k ) are known, find ( Uk+1 , ( Pn ) k+1 , ( Pt ) k+1 ) ∈ IR N × C h n × C h t (( Pn ) k+1 ) such that :            Kk d Uk+1 = Fk + H T n ( Pk+1 ) n + H T t ( Pk+1 ) t Pn k+1 ∈ C h n and ∀P * n ∈ C h n (P * n -( Pn ) k+1 ) T ( Gm n + H n Uk + H n d Uk+1 ) ≥ 0 ( Pt ) k+1 ∈ C h t (( Pn ) k+1 ) and ∀P * t ∈ C h t (( Pn ) k+1 ) (P * t -( Pt ) k+1 ) T (H t Uk + H t d Uk+1 -V obs ) ≥ 0 (53) with Fk = θ 1 F ext + F m Res -1 ∆t M Uk+1 -θ 1 F int (l( Uk )
) and d Uk+1 = Uk+1 -Uk . Kk is the matrix of iterations : Kk = 1 ∆t M + ∆tθ 1 θ 2 K k where K k = ∂F int /∂U is the tangent stiffness matrix which is described in details in [START_REF] Liu | 3d finite element analysis of rubber-like materials at finite strains[END_REF] .

The problem (53) is then formulated in terms of complementarity and, after a condensation on the contact node variables, it is solved by using the Lemke method. This process has been detailled by A.Klarbring and G.Björkmann 11 and by P.Chabrand et al. [START_REF] Chabrand | Various numerical methods for solving unilateral contact problems with friction[END_REF] for a quasistatic problem. The Lemke mathematical programming method is a direct pivoting algorithm which turns out to be very fast 3 .

THE STEADY SLIDING SOLUTION

In the present case, as for disk brakes [START_REF] Moirot | Etude de la stabilité d'un équilibre en présence de frottement de Coulomb, application au crissement des freins à disques[END_REF] , we are interested in instabilities that may occur when a steady frictional regime is desired. Therefore, we propose an algorithm to determine directly that steady sliding solution. The steady sliding problem was mathematically studied for a linear elastic body by P.Rabier et al. [START_REF] Rabier | Existence and local uniqueness of solutions for contact problems with non-linear friction laws[END_REF] (normal compliance law) and by E.B.Pires et al. [START_REF] Pires | The steady sliding problem with nonlocal friction[END_REF] (non-local friction law). Finite dimensional coercive and non-coercive versions of the steady sliding problem with classical unilateral contact conditions were studied by A.Klarbring 9,10 , .

In this section a method is proposed to find a steady sliding equilibrium configuration of a non-linear elastic solid in contact with a flat obstacle that has a constant tangential velocity v obs . The relative tangential velocity of the nodes in contact is equal to -v obs and, consequently, the classical Coulomb friction law simplifies to :

r t v obs = µ|v obs |r n . ( 54
)
For this particular case of frictional contact problems, the direction of the tangential contact stress is known. Thus the problem becomes simpler because the non-smooth character of the friction law is eliminated. We just get the extra relation (54) between the normal and the tangential components of the contact stresses. A new problem is derived from the problem P3, see (48) :

Problem P5 : find (U, P n ) ∈ IR N × C h n such that : F int (U) = F ext + H n P n + µsgn(v obs )H t P n ∀P * n ∈ C h n (P * n -P n ) T (H n U + G) ≥ 0 (55)
where sgn(x) = x/|x| if x = 0. One obtains a sequence of frictionless contact problems by using a fixed point procedure on the frictional tangential contact forces. At each fixed point iteration, the tangential contact forces are updated by using the normal contact forces computed previously and the corresponding problem is :

Problem P6 : (U 0 , P 0 n ), .., (U l-1 , P l-1 n ) being known, find (U l , P l n ) ∈ IR N × C h n such that : F int (U l ) = F ext + H T n P l n + µ sgn(v obs )H T t P l-1 n ∀P * n ∈ C h n (P * n -P l n ) T (H n U l + G) ≥ 0 (56)
where 2 l denotes the value of 2 obtained at the fixed point iteration l. This problem is then linearized similarly to (48). At each iteration of the Newton-Raphson algorithm, the following problem is solved :

Problem P7 : (U l 0 , (P l n ) 0 ), .., (U l k , (P l n ) k ) being known, find (U l k+1 , (P l n ) k+1 ) ∈ IR N × C h n such that : K k dU l k+1 = F ext -F int (U l k ) + H T n (P l n ) k+1 + µ sgn(v obs )H T t P l-1 n ∀P * n ∈ C h n (P * n -(P l n ) k+1 )(H n U l k + H n dU l k+1 + G) ≥ 0 ( 57
)
where 2 l k denotes the value of 2 obtained at the fixed point iteration l and at the Newton-Raphson iteration k and where dU l k+1 = U l k+1 -U l k . The problem P7 (see (57)) is also formulated in terms of complementarity and solved by using the Lemke method.

The initial condition for the fixed point procedure is a steady sliding equilibrium position obtained with µ = 0 or with a smaller value of µ.

The advantage of this method is to get the solutions for various friction coefficients. This is very interesting because the friction coefficient is a key parameter in the occurrence of instabilities. An analysis of the sensitivity of the phenomenon relatively to this parameter will be possible without extra computations.

INSTABILITY OF STEADY SLIDING STATES

In this section we are interested in the occurrence of dynamic instabilities of steady sliding equilibrium states of finite dimensional systems. In other words, we will study the possible existence of growing dynamic solutions in the neighborhood of those equilibrium states. We will show in section 6 that, in the present case, oscillatory (flutter) instabilities are found.

For any given non-vanishing velocity of the flat surface, the sliding velocity of all the nodes in contact at a steady sliding equilibrium configuration is different from zero. We assume, in addition, that the normal reactions of all those contact nodes are strictly positive at the steady sliding state. Thus each contact node is either in (strict) sliding contact (P 0

n i > 0, U 0 n i + G i = 0, P 0 t i = µsgn(v obs )P 0 n i and U 0 t i = 0) or in (strict) separate contact (P 0 t i = P 0 n i = 0 and U 0 n i + G i > 0)
where 2 0 denotes the value of 2 at the steady sliding equilibrium configuration. Then we consider perturbed dynamic solutions in a neighborhood of this steady sliding state, such that the contact state of each contact candidate node is preserved : each node in sliding contact remains in sliding contact (U n i + G i = Un i = 0, P t i = µsgn(P 0 t i )P n i and Ut i = v obs ), and each node not in contact remains so with (P t i = P n i = 0 and U n i + G i > 0). Using the subscript f (free) for the degrees of freedom that correspond to displacement components of node not candidate to contact or not currently in contact, and the subscript n (respectively t) for those that correspond to normal (respectively tangential) displacements of sliding nodes, the dynamical problem may be written :

  M f f M f n M f t M nf M nn M nt M tf M tn M tt     Üf 0 Üt   +   F int f (U f , U t ) F intn (U f , U t ) F intt (U f , U t )   =   0 P n µsP n   ( 58
)
where s = sgn(v obs ). Using the equations of motion corresponding to the normal degrees of freedom of the sliding nodes (n) to eliminate the non vanishing reactions, the equations that govern the smooth dynamic evolution of the structure in that neighborhood are :

M f f M f t M * tf M * tt M * Üf Üt + F int f (U f , U t ) F * intt (U f , U t ) = 0 0 ( 59
)
where F * intt = F intt -µsF intn and M * t2 = M t2 -µsM n2 . If the matrix M * is nonsingular, the second order differential system (59) is regular and the normal displacements of the free contact candidate nodes, as well as the tangential velocities and the normal reactions of the nodes in sliding contact, depend continuously on the initial conditions and on time. Consequently, in a sufficiently small neighborhood of the steady sliding state, the solution of (59) (and ( 58)) for perturbed initial conditions arbitrarily close to U 0 f , U 0 t and U0 f = U0 t = 0 leads to a "smooth" solution of the general system (40), ( 41), (44). The equations of motion (59) are linearized :

M * δ Ü * (t) + K * δU * (t) = 0 ( 60
)
where δU * (t) denotes a small increment of the active displacements, K * a nonsymmetric friction dependent tangent stiffness matrices of dimension N * (≤ N ) constructed as M * .

The main important feature is that if the steady sliding equilibrium state is dynamically unstable for the linearized system (60), it is also unstable for the non-linear system (59) and for the original system (40), (41), and (44). The analysis of the dynamic stability of the system (60) leads classicaly to the generalized eigenproblem :

(λ 2 M * + K * )V * = 0. ( 61
)
The occurrence of a nontrivial eigenvector V * for some λ with a strictly positive real part implies the instability of the steady sliding state: a divergence instability if Im(λ) = 0 (non-oscillatory), a f lutter instability if Im(λ) = 0 (oscillatory).

The generalized eigenproblem (61) is solved by using the Lanczos algorithm [START_REF] Rajakumar | The Lanczos algorithm applied to unsymmetric generalized eigenvalue problem[END_REF] along with the double QR method.

6 SQUEAL FOR A RUBBER/GLASS CONTACT

The waist seal problem

This application adresses the study of the noise generated by the waist seal of the lateral door of a car. This system is fixed to an internal metal sheet and is pressed against the window. The longitudinal and transversal sections of the door are schematically represented in figure 2.

In some cases, this system induces squeal noise. The simulation is conducted on two real waist seals. The first one gives rise to squeal problems and is not used any longer. The second waist seal is a modified version of the first one that has completely removed the noise.

In subsection 6.2, the capabilities of the numerical model are presented on a detailled study of the first waist seal when the glass moves downward. Then, in subsection 6.3, the model is used to compare the behavior of the two waist seals in terms of the possible occurrence of instabilities and of their consequences.

Simulation of a waist seal and stability analysis

The basic parameters for the computation of the waist seal behavior are the geometry, the "loading", the material properties and the coefficient of friction. The geometry is plotted in grey on 3. The mesh is composed of 649 nodes and 538 Q4P0 finite elements.

It has been particulary refined on the contact zone to ensure a minimum of 10 contact nodes. The black triangles cm mark the fixed boundary conditions (only a part of the waist seal is modelled) and the thin rectangle is a schematic representation of the glass window. The loading is applied in two steps : firstly a normal displacement of the obstacle is prescribed and the corresponding equilibrium position is computed for µ = 0. Then the value of the friction coefficient is successively increased, the corresponding steady sliding equilibrium positions being computed for an imposed tangential velocity of the obstacle (80mm/s) either upward or downward. The glass window is supposed to be plane and rectangular and the waist seal to be in a state of plane strain. The Mooney-Rivlin material coefficients are a 1 = 0.293M P a and a 2 = 0.177M P a and the bulk modulus is of 1200M P a. The experimental range of the coefficient of friction is [0.1, 0.8]. This range has been obtained from experiments on small parallelepipedic samples of the waist seal material (with or without varnish) and for various interface conditions (with or without water).

The fixed point algorithm of section 4 is used to compute the steady sliding solutions. Computations are conducted with various friction coefficients (increments of 0.05).

The deformed configuration after the initial normal loading is plotted in black on figure 3. The steady sliding configuration obtained for µ = 0.4 is plotted on figure 4.

For the stability analysis, the eigenvalue spectrum of M * -1 K * is computed for all the steady sliding positions previously obtained. A divergence instability would be characterized by a negative real eigenvalue : it never occurs. Complex eigenvalues correspond to flutter instability. The analysis focuses on the occurrence of instabilities with frequencies belonging to the acoustic range, but some results are given for a larger range. Growth rate (%)

Figure 6: Eigenvalue spectrum of M * -1 K * , µ = 0.4.
The first complex eigenvalue appears for µ = 0.2. The eigenvalue spectrum is plotted in figure 5. The eigenvalue spectrum of M * -1 K * obtained with µ = 0.4 is presented in figure 6.

For µ < 0.285, the frequency of the lowest complex eigenvalue is high (about 44000Hz). On figure 7, the coallescence for µ = 0.285 of two real eigenvalues is shown (frequency about 3250Hz). On the second plot of figure 7, the growth rate of the flutter oscillations is plotted. This complex eigenvalue frequency remains close to 3250Hz for the rest of the range of the friction coefficient.

It can be noted that the frequency associated to the lowest complex eigenvalue is relatively separated from the other ones. The associated eigenmode gives the spatial form of the growing dynamic solution close to the steady sliding position. On figure 8, the deformed configuration is plotted in grey and the deformed configuration plus the eigenmode multiplied by an amplification coefficient is plotted in black. The dynamic behavior is computed by solving the dynamical problem P4 (see ( 53)) taking the steady sliding position and zero velocity as initial conditions. The only initial perturbation is due to the numerical precision in introducing those initial conditions. A slight amount of numerical dissipation is also introduced (θ 1 = θ 2 = 0.55 > 0.50).

The evolution of the global contact forces is plotted on figure 9. During a short time interval, the contact forces seem to remain constant (with P t = µP n ). Then, growing oscillations with the flutter frequency (3250Hz) are perceptible. About 0.25s after the beginning of the dynamical computation, the sliding of some contact nodes becomes intermittent. About 0.03s later, this stick-slip becomes global and is accompanied by global losses of contact (see figure 10). It can be noted on figure 10 that the evolution has reached what looks like a limit cycle. The amplitude of the oscillations of the contact forces is very large and the waist seal hammers the window with a frequency that belongs to the acoustic range. For this computation of the dynamical behavior of the waist seal, the perturbation introduced is very small (only the numerical precision) and a relatively long time is necessary to reach the limit cycle. This time may be reduced by introducing a slightly larger initial perturbation. We have checked that for different forms of this extra small perturbation, the same limit cycle is attained.

For a larger coefficient of friction (µ = 0.8), the growth of the flutter oscillations is larger in the sense that these oscillations become significant for τ = 0.045s (instead of τ = 0.15s for µ = 0.4). The phenomenon is more complex as presented on figure 11. As previously, global losses of contact occur during the oscillations.

An non-linear elastic behavior has been assumed in the theorical model adopted for the waist seal. In reality, this kind of rubberlike material is viscoelastic, and some small numerical damping has been always considered in our computations, by means of the values adopted for θ 1 and θ 2 . It is interesting to check the effect of an increase of the numerical damping on the solution. Using θ 1 = θ 2 = 0.60 instead of θ 1 = θ 2 = 0.55 leads to the results presented on figure 12. The nature of the phenomenon is still the same with an oscillation frequency equal to 3250Hz, but the evolution is more regular and no global losses of contact are observed. This absence of loss of contact has been also observed with µ = 0.4 when θ 1 = θ 2 = 0.60.

We give now some further details of the dynamic evolution in this case (µ = 0.8 and θ 1 = θ 2 = 0.60). We present first the evolution of the tangential velocity of all the contact nodes on a short time interval (τ ∈ [0.25s; 0.2505s]) on figure 13. As expected, a stick-slip phenomenon is observed but it should be noted that it is more complex than a simple succession of global stick and global slip states : the contact nodes do not start to slip simultaneously. We focus now on the the normal and tangential movements of a node located on the edge of the contact zone (node 1 on figure 13). Phase plane plots of that tangential movement are presented in figure 14 (the complete evolution) and in figure 15 (the "limit cycle"). The normal component is presented on figure 16. It can be noted that even for θ 1 = θ 2 = 0.60, the "limit cycle" oscillation has a short period of loss of contact; for example, one cycle of the node 1 involves the following sequence of states : stick, slip, loss of contact, slip and back to stick. The results obtained with the window moving upward are very similar to the previous one, but with a lower flutter frequency (about 1980Hz).

Comparison of the analysis for two different waist seals

Because of the lack of detailed experimental results, the only way to give a rough validation of the model has been to compare the results of the stability analyses of two waist seals, one of which is known to be noisy and the other one to be "silent". Those two waist seals are flocked, and with this flock the friction coefficient is experimentally found to belong to [0.15, 0.45]. The waist seal simulated in the previous subsection was known to be noisy. We will conduct now the analysis of the modified "silent" geometry. This second geometry is shown to be less prone to flutter instabilities than the first one. The new geometry is given on figure 17. The eigenvalue spectra of M * -1 K * for the upward movement of the window with µ = 0.4 are presented in figure 19 (first geometry) and in figure 20 (second geometry). The eigenmode associated to the lowest flutter frequency is plotted on figure 18. For the first case, the frequency corresponding to the complex eigenvalue is of 1950Hz with a growth rate of 3.5% and for the second case, it is of 18890Hz with a growth rate of 0.06%. For the second geometry, squeal noise is not expected, because, for the computed high frequency of flutter (which in any case would not be in the acoustic range) and for the corresponding Growth rate (%) small growth rate, one may expect that any small amount of material dissipation will cancel out the instability. For the window moving downward, the eigenvalue spectrum of M * -1 K * , corresponding to the first case has been given on figure 6. The one corresponding to the second waist seal is given on figure 21. The frequencies corresponding to the first complex eigenvalue are similar in both cases (3250Hz for the first geometry, 3050Hz for the second one) but the growth rate is six times smaller for the second geometry (0.35% instead of 2.4%).

The second case involves also two other frequencies in the acoustic range (5240Hz and 6810Hz) but their growth rates are also smaller than the one of the first geometry. To check the possible consequences of flutter in this second case, dynamic analyses were performed with the numerical parameters used earlier to obtain figure 22 (namely θ 1 = θ 2 = 0.55). The result given on figure 22 shows that no oscillating solution emerges.

The first geometry is more prone to flutter instability then the second one and, using the same material and numerical parameters, important oscillations are numerically obtained for the first "noisy" geometry whereas no oscillations occur in the numerical tests with the second "silent" geometry.

For sufficiently large friction coefficient (for µ = 0.8) the second geometry may also lead to squeal (the first one was unstable for both µ = 0.4 and µ = 0.8). The lowest flutter frequency obtained with µ = 0.8 is of 2965Hz and the corresponding growth rate is of 3.35%. The dynamic solution plotted on figure 23 should be compared with the one obtained for the first geometry and given on figure 11.

Figure 2 :

 2 Figure 2: Schematic description of the lateral door.

Figure 3 :

 3 Figure 3: Initial normal loading, wiper compression with µ = 0.

Figure 4 :Figure 5 :

 45 Figure 4: Glass downward: steady sliding position for µ = 0.4.

Figure 7 :

 7 Figure 7: Coupling of two real eigenvalues, µ ∈ [0.25 : 0.3].

Figure 8 :

 8 Figure 8: Eigenmode associated to the first flutter frequency, µ = 0.4, first version.

Figure 9 :Figure 10 :

 910 Figure 9: Evolution of the contact forces for µ = 0.4.

Figure 11 :Figure 12 :

 1112 Figure 11: Evolution of the contact forces for µ = 0.8, θ 1 = θ 2 = 0.55.

Figure 13 :Figure 14 :

 1314 Figure 13: Evolution of the tangential velocity of the contact nodes with τ ∈ [0.25s; 0.2505s], the window velocity is of 80mm/s.

Figure 15 :Figure 16 :

 1516 Figure 15: Limit cycle, τ ∈ [0.2; 0.25].

Figure 17 :

 17 Figure 17: Steady sliding position, window downward, µ = 0.4, second waist seal.

Figure 18 :

 18 Figure 18: Eigenmode associated with the lowest flutter frequency, µ = 0.4, second waist seal.

Figure 19 :Figure 20 :

 1920 Figure 19: Eigenvalues spectrum of M * -1 K * , µ = 0.4, window moving upward, first waist seal.

Figure 21 :

 21 Figure 21: Eigenvalues spectrum of M * -1 K * , µ = 0.4, window moving downward, second waist seal.

6. 3 . 2

 32 Window moving downward (µ = 0.4 and θ 1 = θ 2 = 0.55)

Figure 22 :Figure 23 :

 2223 Figure 22: Evolution of the contact forces for µ = 0.4, θ 1 = θ 2 = 0.55, second waist seal.

  n i and H t i are two vectors of dimension N . The matrix H n (respectively H t ) of dimension N c × N puts together the N c vectors H n i (respectively H t i ). We denote also by H the matrix of dimension 2 * N c × N such that H T = H T n , H T t . The sets C n and C

t (R n ) are discretized in C h n and C h t (P n ) :

ACKNOWLEDGEMENT

Research supported by Renault in the project "Numerical study of the friction between a glass window and a wiper".


Clément Acquitter

Benjamin Lemasson

Elizabeth Moyal

Soléakhéna Ken

IMPACT OF MULTICENTER HARMONIZATION ON THE DETECTION OF RADIONECROSIS IN MULTIMODAL MRI RADIOMICS STUDY

Keywords: radiomics, multi-modal MRI, harmonization, glioblastoma, radionecrosis

whether they are published or not. The documents may come

INTRODUCTION

Hypo-fractionated stereotactic radiotherapy (hFSRT) is an effective treatment for recurrent GBM. Because a high gradient dose is delivered, radionecrosis may occur and complicate MRI interpretation. Radiomics analysis can be used to detect radionecrosis-related features [1]. However, the multi-centric nature of the data may affect the performance of such models. The impact of feature harmonization on the detection of radionecrosis is detailed in this study.

MATERIAL AND METHODS

Multi-center clinical MRI database includes 75 examinations of T1w, T1w post-contrast, T2w, and FLAIR images from 6 centers. Brain tumor segmentation was performed thanks to a deep learningbased algorithm implemented in the BraTs Toolkit [START_REF] Kofler | Frontiers in Neuroscience[END_REF]. A manual correction was applied to retrieve a 3-label segmentation corresponding to the following areas: hyper-enhanced tumor, peritumoral edema, resection cavity (Fig. 1). A total of 1200 radiomics features were extracted from each examination (100 first-order, textural, and shape features for each label and each MR modality). The predictive ability of 18 radiomics models (6 features selection, namely chi2, F-classification and Mutual Information with respectively 5 and 10 features and 3 classification methods, namely logistic regression (LR), SVM and Random Forest (RF)) were evaluated with and without features harmonization with ComBat. [START_REF] Fortin | [END_REF].

RESULTS

For the LR, SVM and RF classification, mean area under the curve (AUC) scores was respectively 0.67 ± 0.04, 0.51 ± 0.17 and 0.70 ± 0.14 without harmonization versus 0.75 ± 0.05, 0.75 ± 0.03 and 0.77 ± 0.03 with harmonization. Heatmaps show increased performance for most of the models trained on the harmonized features (Fig. 2). The best score increase was obtained for the SVM classification method (+105%). Models based on RF showed slight performance decreases (from -0.4% up to -15%) with the harmonization process. Overall, the harmonization of radiomics features increases the precision of the detection of radionecrosis while decreasing the variability induced by the selection or the classification choices (0.75 ± 0.04 vs. 0.67 ± 15; mean ± SD).

CONCLUSION

Our results show that multicenter data harmonization increases the performance of radiomics models dedicated to the detection of radionecrosis in patients with recurrent GBM treated by hFSRT.

COMPLIANCE WITH ETHICAL STANDARDS

Patients gave their consent in agreement with the Ethic committee.

ACKNOWLEDGMENTS

This work is financially supported by the Project SIGN'IT.

Fig. 1 .

 1 Fig. 1. Labels of hyper-enhanced tumor (yellow), peritumoral oedema (green) and resection cavity (pink) are mapped over T1w (A), FLAIR (B), T1w post-contrast (C) and T2w (D) MRI.

Fig. 2 .

 2 Fig. 2. Classification scores (AUC) for each radiomics model without (left) and with harmonization (right)