Evang Éline Capiez-Lernout
email: evangeline.capiez-lernout@u-pem.fr

Christian Soize
email: christian.soize@u-pem.fr

Moustapha Mbaye
email: moustapha.mbaye@turbomeca.fr

COMPUTATIONAL GEOMETRICALLY NONLINEAR VIBRATION ANALYSIS OF UNCERTAIN MISTUNED BLADED DISKS

The recent improvements in turbomachinery design requires the analysis of exceptional operating regime of bladed disks corresponding to geometrical nonlinear effects induced by the large displacements/deformations. In addition, the random nature of the mistuning has also to be modeled. First, a mean nonlinear reduced-order model of the tuned bladed disk is explicitly constructed in the context of the finite element method. The investigation is then devoted to the modeling of the mistuning through the nonparametric probabilistic approach extended to the nonlinear geometric context. The stochastic nonlinear equations are solved in the time domain using the Monte Carlo numerical simulation coupled with advanced arc-length methods adapted to high nonlinear response levels. Finally, the methodology is applied through a numerical example of a bladed disk and a nonlinear analysis is performed in both time and frequency domain. u Displacement vector (mean computational model) ũ Observation vector (mean NL-ROM) [A] Finite element matrices (Roman capital letters) [Φ] Modal matrix

NOMENCLATURE

Frequency band of analysis s Frequency band of excitation Ê Real vector space Spatial discretization of the load f External load vector f NL (u) Nonlinear restoring force q vector of the generalized coordinates (mean NL-ROM)

INTRODUCTION

There are many aspects to be taken into consideration when dealing with turbomachinery bladed disks. The predictive computational models are more and more sophisticated, including the modeling of complex phenomena occurring in such structures. It is well known, that the manufacturing tolerances, the dispersion of materials, which arouses in bladed disks structures create uncertainties on the geometry, on the boundary conditions and on the material properties, which breaks the natural cyclic symmetry of such structures. The main consequence results in strong vibrations combined to spatial localization in the dynamic forced response of the blades [START_REF] Wei | Localization phenomena in mistuned assemblies with cyclic symmetry part ii: Forced vibrations[END_REF]. The random character of such phenomenon referred to as mistuning requires the use of adapted probabilistic model of uncertainties in the computational models [START_REF] Bladh | Component-mode-based reduced order modeling techniques for mistuned bladed disks -part i: Theoretical models[END_REF][START_REF] Rivas-Guerra | Identification of mistuning characteristics of bladed disks from free response data -part ii[END_REF][START_REF] Mignolet | Nonparametric stochastic modeling of linear systems with prescribed vari-ance of several natural frequencies[END_REF][START_REF] Capiez-Lernout | Nonparametric modeling of random uncertainties for dynamic response of mistuned bladed disks[END_REF]. Another essential aspect to be taken into account concerns the modeling of the different nonlinearity sources such as nonlinear contact interfaces [START_REF] Petrov | Analysis of sensitivity and robustness of forced response for nonlinear dynamic structures[END_REF][START_REF] Petrov | Analysis of flutter-induced limit cycle oscillations in gas-turbine structures with friction, gap, and other nonlinear contact interfaces[END_REF][START_REF] Laxalde | Qualitative analysis of forced response of blisks with friction ring dampers[END_REF] or geometrically nonlinear effects [START_REF] Vakakis | Dynamics of a nonlinear periodic structure with cyclic symmetry[END_REF][START_REF] Sarrouy | Global and bifurcation analysis of a structure with cyclic symmetry[END_REF][START_REF] Grolet | Vibration analysis of a nonlinear system with cyclic symmetry[END_REF][START_REF] Grolet | Free and forced vibration analysis of a nonlinear system with cyclic symmetry: Application to a simplified model[END_REF]. At the same time, the progression of the computational capabilities, which includes the possibility of using parallel computations, has been a motivation to develop computational strategies adapted to large finite element model of industrial bladed disks. Various researches concerning reduced-order modeling techniques adapted to the mistuned linear analysis of bladed disks [START_REF] Bladh | Component-mode-based reduced order modeling techniques for mistuned bladed disks -part i: Theoretical models[END_REF][START_REF] Yang | A reduced-order model of mistuning using a subset of nominal modes[END_REF] have been used for robust design optimization purpose [START_REF] Castanier | Using intentional mistuning in the design of turbomachinery rotors[END_REF][START_REF] Choi | Optimization of intentional mistuning patterns for the reduction of the forced response effects of unintentional mistuning: Formulation and assessment[END_REF][START_REF] Mbaye | Robust analysis of design in vibration of turbomachines[END_REF]. Adapted reduced-order techniques have also been proposed in [START_REF] Krack | Reduced order modeling based on complex nonlinear modal analysis and its application to bladed disks with shroud contact[END_REF][START_REF] Petrov | A high-accuracy model reduction for analysis of nonlinear vibrations in structures with contact interfaces[END_REF] for the nonlinear dynamical analysis of turbomachines.

In the present case, we are interested in the geometrically nonlinear analysis of mistuned bladed disks. The recent improvements in turbomachinery design requires the analysis of exceptional operating regime of bladed disks for which large displacements/deformations can occur. The case of severe loading is investigated in the context of elasto-dynamics. Such situation is equivalent to nearly unstable cases induced by aerodynamic coupling yielding flutter and thus very low damping levels. In such case, the linearized elasto-dynamic theory can not be used anymore because the geometrically nonlinear effects induced by the large deformations and the large displacements are very strong and need to be taken into account in the modeling.

The methodology is developed in the context of large finite element computational model of bladed disks, requiring the construction of adapted nonlinear reduced-order computational models (NL-ROM). The developments of such NL-ROM requires the selection of an appropriate deterministic basis for the representation of the nonlinear dynamic response. More generally, this projection basis can be obtained by using the proper-orthogonal decomposition method or by an eigenvalue analysis. A mixed formulation combining these two approaches has been proposed in [START_REF] Capiez-Lernout | nonlinear stochastic dynamics for post-buckling analysis of uncertain cylindrical shells[END_REF]. Strategies based on the selective choice of the most appropriate basis are reviewed in [START_REF] Mignolet | A review of indirect/non-intrusive reduced order modeling of nonlinear geometric structures[END_REF]. The parameters of the NL-ROM of the tuned structure can then be either deduced using the STEP procedure (which is based on the smart non-intrusive use of standard commercial finite element codes) [START_REF] Mignolet | A review of indirect/non-intrusive reduced order modeling of nonlinear geometric structures[END_REF][START_REF] Muryavov | Determination of nonlinear stiffness with application to random vibration of geometrically nonlinear structures[END_REF][START_REF] Mignolet | Stochastic reduced order models for uncertain geometrically nonlinear dynamical systems[END_REF] or from explicit construction as shown in [START_REF] Capiez-Lernout | Computational stochastic statics of an uncertain curved structure with geometrical nonlinearity in threedimensional elasticity[END_REF] in the context of three-dimensional solid finite elements. Once the NL-ROM of the tuned bladed disk is established, uncertainties can be implemented through the nonparametric framework [START_REF] Soize | A nonparametric model of random uncertainties for reduced matrix models in structural dynamics[END_REF][START_REF] Soize | Stochastic Models of Uncertainties in Computational Mechanics[END_REF], which has been extended in the nonlinear geometric case [START_REF] Mignolet | Stochastic reduced order models for uncertain geometrically nonlinear dynamical systems[END_REF]. Such probabilistic approach is able to simultaneously capture both system-parameter uncertainties and model uncertainties. It relies on the dedicated construction of a probability model related to the tuned NL-ROM operators, using the maximum entropy principle (MaxEnt).

The paper is organized as follows. Section 1 gives a complete description of the computational methodology and of its numerical aspects allowing the nonlinear dynamic analysis of mistuned bladed-diks to be performed. The mean NL-ROM corresponding to the tuned structure is explicitly constructed following the approach presented in [START_REF] Capiez-Lernout | Computational stochastic statics of an uncertain curved structure with geometrical nonlinearity in threedimensional elasticity[END_REF]. Once the mistuning is modeled through the nonparametric probabilistic approach, the numerical strategy concerning the resolution of the set of nonlinear coupled differential equations is discussed. Section 2 is devoted to a numerical example, for which the geometrically nonlinear effects are analyzed and quantified in both tuned and mistuned cases.

METHODOLOGY

This Section is devoted to the construction of a methodology for the nonlinear mistuning analysis occurring in rotating bladed disks structures. In the present research, the bladed disks under consideration are assumed (1) to be made up of a linear elastic material and (2) to undergo large displacements and large deformations inducing geometrical nonlinearities.

Nonlinear dynamics of a tuned bladed disk

The tuned bladed disk structure has an N-order cyclic symmetry. Thus, the geometrical domain, the material constitutive equations and the boundary conditions related to the generating sector are invariant under the 2π N rotation around its axis of symmetry. Moreover, the bladed disk undergoes a rotational motion around the axis of symmetry with constant angular speed Ω. By choosing a total Lagrangian formulation, the dynamical equations are expressed in the rotating frame of an equilibrium configuration considered as a prestressed static configuration.

The mean (or nominal) computational model of the tuned bladed disk, which is constructed by the Finite Element Method (FEM) is written as [START_REF] Desceliers | Nonlinear viscoelastodynamic equations of three-dimensional rotating structures in finite displacement and finite element discretization[END_REF]:

[M ] ü + ([C(Ω)] + [D ]) u + [K 1 (Ω)]u + f NL (u) = f , [K 1 (Ω)] = [K e ] + [K g ] + [K c (Ω)] , (1)
in which the Ê n -vector f is the external load representing for in- stance the unsteady pressures applied to the blades, and where the Ê n -vector u corresponds to the finite element discretization of the unknown displacement field. In Eq. ( 1), the matrices [M ],

[D ], [K g ] and [K e ] are the mass, damping, geometrical stiffness and elastic stiffness real matrices with positive definiteness property. The rotational effects are taken into account through the gyroscopic coupling matrix [C(Ω)] and the centrifugal stiffness matrix [K c (Ω)], which have respectively antisymmetry property and negative definiteness property. It should be noted that all these matrices are also N-block circulant matrices [START_REF] Davis | Circulant matrices[END_REF] since the structure has an N-order cyclic symmetry. The geometrical nonlinearities effects are taken into account through the Ê n -vector f NL (u) which includes the exact and complete quadratic and cubic terms given by the three-dimensional nonlinear geometric elasticity. Furthermore, the centrifugal effects are assumed to be sufficiently small so that the linear stiffness matrix [K 1 (Ω)] is positive definite, yielding only stable dynamical systems to be considered.

In the present case, the presence of the geometric nonlinearity naturally yields the nonlinear equations to be solved in the time domain, the frequency content of the nonlinear response being a posteriori analyzed by Fourier transform. Concerning the external load, a usual harmonic excitation would be inappropriate because the nonlinear equations should be solved for each harmonic excitation considered. The external load is then defined in the time domain corresponding to a uniform sweep of a chosen frequency band of excitation. Let s =s ∪ s be the frequency band of excitation with central frequency s∆ ω and bandwidth ∆ω defined by s = [(s-1/2)∆ω, (s+ 1/2)∆ω]. The external load is written as

f(t) = f 0 g(t) , (2)
in which f 0 is a coefficient characterizing the global load intensity, and where is an Ê n -vector corresponding to the spatial discretization of the load. In Eq. ( 2), the function g(t) is cho-

sen as g(t) = ∆ω π sinc π ( t∆ω 2π ) cos(s ∆ω t), where x → sinc π (x)
is the special function defined by sinc π (x) = sin(π x)/(π x). It should be noted that all the frequencies of the frequency band of excitation are simultaneously excited so that only one nonlinear time-domain analysis is carried out. With such time-domain excitation, a forced-response problem is considered and not a time evolution problem with initial conditions. The considered forced-response problem is thus approximated by an equivalent time-evolution problem with zero initial conditions over a finite time interval, which includes almost all of the signal energy of the excitation. The use of the cyclic symmetry property by decomposing the nonlinear response according to its harmonic components is not considered in the present case because all the harmonic components are coupled through the geometric nonlinearity and does not allow the problem on a single rotor sector to be solved. Moreover, the full computation of the nonlinear solution of Eq. ( 1) can not be reasonably performed when dealing with realistic models of bladed disks corresponding to a large number of DOF. More specifically, large dimensional computational models require to develop reduced-order model strategies adapted to this geometric nonlinear context (see [START_REF] Capiez-Lernout | nonlinear stochastic dynamics for post-buckling analysis of uncertain cylindrical shells[END_REF][START_REF] Muryavov | Determination of nonlinear stiffness with application to random vibration of geometrically nonlinear structures[END_REF][START_REF] Capiez-Lernout | Computational stochastic statics of an uncertain curved structure with geometrical nonlinearity in threedimensional elasticity[END_REF] and [START_REF] Mignolet | A review of indirect/non-intrusive reduced order modeling of nonlinear geometric structures[END_REF] for a complete overview). Let be a given vector basis represented by the (n × P) real matrix [Φ]. The nonlinear response u is expanded as

u = [Φ] q , (3)
in which q is the Ê P -vector of the generalized coordinates. Re- placing Eq. ( 3) in Eq. ( 1) yields a nonlinear reduced set of P coupled differential equations for which all linear, quadratic and cubic terms have to be known. In the present research, the construction of the operators of such mean NL-ROM is explicitly carried out in the context of the three-dimensional finite element method. The finite elements are isoparametric solid finite elements with 8 nodes using a numerical integration with 8 Gauss integration points. The elementary internal forces projected on the chosen vector basis are numerically constructed for each finite element before performing its assembly and computing all linear and nonlinear reduced operators. The detailed procedure, which also uses the symmetry properties of the linear and nonlinear reduced-operators combined with distributed computations, can be found in [START_REF] Capiez-Lernout | Computational stochastic statics of an uncertain curved structure with geometrical nonlinearity in threedimensional elasticity[END_REF]. It should be noted that each type of reduced operator is separately modeled, keeping open the possibility of implementing uncertainties issued from independent physical sources.

Nonlinear dynamics of a mistuned bladed disk

The random nature of the mistuning is then considered by implementing the nonparametric probabilistic approach, which presents the ability to include both system-parameter uncertainties and model uncertainties (see [START_REF] Soize | Stochastic Models of Uncertainties in Computational Mechanics[END_REF] for a complete review on the subject). It consists in replacing the operators of the mean NL-ROM by random operators, whose probability distribution is derived from the maximum entropy principle. In the present case, the mass, damping, geometrical stiffness and the opposite of the centrifugal stiffness are positive definite operators whose probability model is constructed from the usual nonparametric probabilistic approach [START_REF] Soize | A nonparametric model of random uncertainties for reduced matrix models in structural dynamics[END_REF]. The linear elastic, quadratic and cubic stiffness reduced operators are reshaped in a positive-definite operator as shown in [START_REF] Mignolet | Stochastic reduced order models for uncertain geometrically nonlinear dynamical systems[END_REF], which allows the nonparametric probabilistic theory to be extended to the geometric nonlinear case. The stochastic model of the gyroscopic coupling is constructed with a similar method as proposed in [START_REF] Mbaye | A reducedorder model of detuned cyclic dynamical systems with geometric modifications using a basis of cyclic modes[END_REF] for the aero-elastic coupling. The dispersion of each random operator is then characterized by one scalar parameter. Consequently, the mistuning level of the bladed disk is entirely controlled by the Ê 6 -vector

= (δ M , δ D , δ C , δ K g , δ K c , δ K ).

Numerical computations

The solution of the stochastic NL-ROM is calculated using the Monte Carlo numerical simulation. For each realization, a set of P nonlinear coupled differential equations is considered and solved with the Newmark method, for which the averaging acceleration scheme, known to be unconditionally stable is used. With this solver, a set of nonlinear equations whose solution is denoted by the Ê P -vector Õ has to be solved at each sampling time. Such computation is mainly addressed by the fixed point method because the iterative scheme does not require the evaluation of the tangential stiffness matrix. Nevertheless, when the algorithm does not converge, it is replaced by the Crisfield arc-length method [START_REF] Crisfield | Non-linear finite element analysis of solids and structures[END_REF]. Such algorithm introduces an additional scalar unknown that multiplies the right-hand side of the nonlinear equation. It is solved step by step, each incremental step being characterized by a given arc-length. For a given step, an iterative scheme requiring one evaluation of the tangential stiffness matrix allows a solution (Õ, ) to be computed. An adaptive arc-length, depending on the number of iterations necessary to obtain the convergence of the preceding increment is also implemented. Furthermore, since parameter has to be controlled to be equal to 1, the state of the algorithm corresponding to the preceding increment has to be stored. Even if the procedure is time consuming, its main advantage concerns its capability of capturing high non-linear mechanical behaviors.

NUMERICAL EXAMPLE Description of the mean computational model

The structure under consideration is a bladed disk constituted of a disk and 24 blades. The disk is made of a homogeneous and isotropic material with constant thickness 0.005m, inner radius 0.035m, outer radius 0.1m, mass density 7860Kg.m -3 , Poisson ratio 0.25 and Young modulus 2 × 10 11 N.m -2 . Fixed conditions are applied along the internal boundary defined by the inner radius. Each blade is made of the same homogeneous and isotropic material as the disk one with length 0.07 m, width 0.0085m, linear decreasing thickness from 0.005 m to 0.001 m. The structure is in rotation around its revolution axis with a constant angular speed Ω = 30 000 rpm. Since the dynamic analysis is carried out in the rotating frame of the structure, the rigidbody motion due to the rotation of the structure corresponds to a fixed boundary condition at the inner radius of the structure. A damping model is added for the bladed disk, corresponding to a hysteretic model with a mean loss factor 2 × 10 -7 , which corresponds to a modal damping around 0.001 for the low frequency band of interest defined by = [0 , 3 900] Hz.

The full finite element model, shown in Fig. 1, is built with 288 8-nodes solid finite elements and is thus constituted of 2 592 degrees of freedom (see Table 1 for the detailed data related to the mesh of the structure).

Nonlinear analysis of the tuned bladed disk

Tuned linear eigenfrequencies The natural eigenfrequencies of the linear tuned bladed disk are calculated first by using the cyclic symmetry property of the structure [START_REF] Thomas | Dynamics of rotationnally periodic structures[END_REF][START_REF] Bossavit | Boundary value problems with symmetry and their approximation by finite elements[END_REF]. Figure 2 displays these eigenfrequencies as a function of its circumferential wave number.

Description of the external load

Concerning the modeling of the external load according to Eq.( 2), the frequency band of excitation is chosen with regard to Fig. 2 such that s = [3000 , 3207] Hz, which corresponds to s = 15 and ∆ω = 1 300 rad.s -1 . The spatial distribution of the load is only concerned with point excitations located at the end of each blade along the transverse direction. Vector characterizing the spa-

FIGURE 2. NATURAL FREQUENCIES WITH RESPECT TO CIR-CUMFERENTIAL WAVE NUMBER

tial distribution from blade to blade has then 24 non zeros components, which are uniformly distributed from blade to blade. In that case when the tuned linear case is considered, there are contributions on all the circumferential wave number. The external load is thus considered as a function of f 0 . Concerning the numerical sampling, the initial instant of integration t ini and the total time duration T are chosen according to [START_REF] Soize | Medium frequency linear vibrations of anisotropic elastic structures[END_REF] such that t ini = -0.029 s (corresponding to a zero value of function g(t)) and T = 0.36 s. Given combined with the Shannon theorem yields the number n t of time steps to be n t = 2 800. The frequency resolution is then δ ν = 2.78 Hz. Let ĝ(ω) be the Fourier transform of function g(t). Figure 3 shows the graphs t → g(t)

and ν → ĝ(2πν), in which ν = ω 2π .

Convergence analysis of the mean NL-ROM Concerning the choice of the vector basis for the construction of the mean NL-ROM, a strategy would consist in using the proper orthogonal decomposition method (POD). This would require the knowledge of the full reference solution, which is practically not possible to get, when considering realistic bladed disk models involving large number of degrees of freedom (DOF). In the present case, the nonlinear equations are solved in the subspace spanned by the usual linear basis constituted of the P modal shapes related to the first increasing natural eigenfrequencies, according to [START_REF] Capiez-Lernout | Computational stochastic statics of an uncertain curved structure with geometrical nonlinearity in threedimensional elasticity[END_REF]. Zero initial conditions are used. The observation is given by the Ê N -vector ũP corresponding to the out-plane dof located at the end of each blade. An incremental convergence analysis is carried out for a maximum load of f 0 || || = 3.512 N, corresponding to an high rate of geometrical nonlinearity. Such high loading can be interpreted as equivalent to a situation for which the damping of the bladed disk structure would reach very small values. Such extreme situations are realistic when flutter occurs and yields nearly unstable situations. Let Conv(P) be the function defined by

Conv(P) = t ini +T t ini || ũP (t)|| 2 dt . ( 4
)
Figure 4 displays the graph P → Conv(P)/Conv(150) in a logarithmic scale. It can be seen that P = 90 yields a reasonable convergence. From now on, the calculations will be carried out for P = 100.

Nonlinear time domain analysis From now on, the converged solution corresponding to the observation issued from the mean NL-ROM is denoted by ũ(t). For clarity, when confusion is possible, superscripts LIN and NONLIN will be added for distinguishing the linear case from the geometric nonlinear one. Let be the load rate such that r = 100% when intensity load is

f 0 || || = 3.512 N.
The intensity loads considered are then described by r f 0 || || with r ∈]0 , 1].

Being interested in the blade yielding the highest vibration amplitude, let j 0 = arg max j max t ũNONLIN j (t) calculated for a load rate r = 100%. The observation v(t) corresponding to the selected blade out-plane displacement is defined by v(t) = ũ j 0 (t). Figures 5 and6 display the graphs t → v LIN (t) (upper graph) and t → v NONLIN (t) (lower graph) for load rates r = 20% and r = 100% corresponding to significant levels and very high levels of geometric nonlinear effects. On these two figures, it is seen that the geometric nonlinearities induce a blade stiffening characterized by a reduction of the vibration amplitudes of the blades with respect to the linear case. This stiffening is also combined with a strong irregularity of the blade response shape over time, which shows an enrichment of the frequency content, which has to be quantified. larges with increasing load rate. It should be noted that there also exist higher frequencies excited through this nonlinearity which are not observed in the chosen band of analysis . As observed on the graphs, the amplitude levels from linear and high nonlinear cases are drastically different.

Nonlinear frequency domain analysis

Nonlinear analysis of the mistuned bladed disk Mistuning implementation

The mean NL-ROM used above is re-used for considering the mistuned bladed disk. The stochastic NL-ROM is defined from the mean NL-ROM by replacing the operators of the mean NL-ROM by random operators according to the nonparametric probabilistic theory [START_REF] Soize | Stochastic Models of Uncertainties in Computational Mechanics[END_REF]. Since the NL-ROM is constructed by modal analysis without substructuring techniques, the uncertainties are not considered to be independent from one blade to another one, which restricts the analysis to the class of integrated bladed disks, manufactured from a unique solid piece of material. In order to simplify the analysis, the extension being straightforward, only the linear elastic part is considered to be uncertain. In that case, the reduced positive-definite operator

[Φ] T [K e ] [Φ], which can be writ- ten such that [Φ] T [K e ] [Φ] = [L] T [L] is replaced by the random operator [L] T [G(δ )] [L], in which [G(δ )
] is a random operator whose probability distribution is entirely described in [START_REF] Soize | Stochastic Models of Uncertainties in Computational Mechanics[END_REF]. The dispersion parameter δ is then a scalar parameter. According to the nonparametric probabilistic theory, it allows the uncertainty level of the stochastic NL-ROM to be controlled. In the present case, numerical simulations are carried out for dispersion level δ = 0.02 (small mistuning) and δ = 0.1 (moderate mistuning) and with r = 100%.

Nonlinear time domain analysis

The nonlinear mistuning analysis is first considered in the time domain. For fixed t, observation v(t) is replaced by the random variable V (t). The numerical simulations are carried out with n s = 1024 Monte Carlo realizations, which is sufficient to get the convergence of {V 2 (t)}, in which is the mathematical expectation. Let V ∞ be the random variable defined by Figure10 displays the graph of the probability distribution function v → p V NONLIN ∞ (v), which is computed using the kernel density estimation [START_REF] Bowman | Smoothing Techniques for Data Analysis[END_REF]. It has to be compared with the value v = 3.665 × 10 -4 m related to the nonlinear analysis of the tuned bladed disk (see Fig. 6). When δ increases from 0.02 to 0.1, the coefficient of variation of the observation is doubled from 0.027 to 0.059, yielding a broader range taken by the Monte Carlo realizations. Moreover, a significant loss of symmetry of the probability distribution function, quantified by a variation of the skewness from -0.07 to 0.86 can be seen. It is characterized by a shift of the support to higher values, yielding a few realizations of the observation to be increased of 40% with respect to the nonlinear tuned case.

V ∞ = max t∈[t ini ;t ini +T ] V (t) . (6)
Nonlinear frequency domain analysis A Fourier transform is then conducted in order to perform the nonlinear analysis in the frequency domain. For fixed ω, let Y (ω) be the random dynamic magnification factor defined by

Y (ω) = W (ω) max ω∈ w(ω) , (7)
in which W (ω) is the random variable similar to observation w(ω) in the nonlinear tuned case. It should be noted that Y (ω) is a normalized random quantity and that the amplitude levels between the geometric nonlinear case and the linear case are about an order smaller. Figures 11 and12 display the graphs of the confidence region of ν → Y LIN (2πν) calculated with a 95% probability level for r = 100%, δ = 0.01 and δ = 0.2. Figures 13 and14 the similar graphs for random observation ν → Y NONLIN (2πν). Comparing the linear case with the nonlinear case shows that the effects of uncertainties increase with the geometric nonlinearity. For fixed δ , it can be seen on these figures that the confidence region related to the peak located around s is very sensitive to uncertainties, when considering the geometric nonlinear case, yielding a broad confidence region and some realizations with 40% (or 30%) amplification effect with δ = 0.02 (or δ = 0.1). A small softening effect is also observed with a 0.5% decreasing shift of the frequency resonance. The uncertainties are also propagated on the secondary peaks excited through the geomet-

CONCLUSION

The paper has presented an analysis of the geometrical nonlinear effects of uncertain mistuned bladed disk, corresponding to the particular case for which aerodynamic coupling yields flutter, that is to say nearly unstable situations with very low damping, which is represented by an high loading level. Firstly, a nonlinear dynamic analysis of the tuned bladed disk is proposed through the construction of a NL-ROM. The effects of the geometric nonlinearities, corresponding to a quick stiffening of the structure when reaching a critical level of loading, are quantified in both time domain and frequency domain. The dynamical response of the blades is also investigated outside the frequency domain of excitation. The linear response has its energy concentrated in the frequency domain of excitation whereas the nonlinear response is spread over a larger frequency domain. Then, the nonlinear analysis of the mistuned structure is proposed in the context of integrally bladed disks, assuming the uncertainties from one blade to another one to be dependent. Compared to the linear mistuned case, the nonlinear mistuned response predictions yield low vibration amplitudes. Nevertheless, for a given mistuning rate, the nonlinear mistuned response predictions display broader confidence regions and are much more sensitive to uncertainties. Such decrease of robustness with respect to uncertainties propagates quickly in the whole frequency band of analysis. This yields a complex vibratory situation, giving rise to frequency ranges for which secondary resonances appear with large confidence regions.

ACKNOWLEDGMENT

This work was supported by the DGA (French defence procurement agency) in the context of the TURBODYNA project (project number ANR-13-ASTR-0008-01) related to the ANR ASTRID research program (specific support scheme for research works and innovation defence).

  g(t) Time domain of the load ĝ(ω) Frequency domain of the load j 0 , k 0 Blade number n number of degrees of freedom n t Number of time steps r Load rate s∆ω Centered frequency of s t ini Initial instant of integration v(t) Time domain observation (mean NL-ROM) w(ω) Frequency domain observation (mean NL-ROM) N Order of the cyclic symmetry P Dimension of the NL-ROM T Total duration time V (t) Time domain observation (stochastic NL-ROM) V ∞ Time domain observation (stochastic NL-ROM) W (ω) Frequency domain observation (stochastic NL-ROM) Y (ω) Magnification factor (stochastic NL-ROM) δ ν Frequency resolution δ Dispersion parameter ν Frequency (Hz) ω min Lower bound of s ω Pulsation (rad.s -1 ) ∆ω Width of s Ω Angular speed (rotational motion)

FIGURE 1 .

 1 FIGURE 1. FINITE ELEMENT MESH.

2 FIGURE 3 .

 23 FIGURE 3. REPRESENTATION OF THE EXTERNAL LOAD IN THE TIME AND FREQUENCY DOMAIN: GRAPH OF t → g(t) (UP-PER GRAPH) and ν → ĝ(2πν) (LOWER GRAPH).

FIGURE 4 .FIGURE 5 .

 45 FIGURE 4. CONVERGENCE ANALYSIS OF THE NL-ROM OF THE TUNED BLADED DISK: GRAPH OF FUNCTION P → Conv(P)

FIGURE 6 .FIGURE 7 .FIGURE 8 .

 678 FIGURE 6. TIME DOMAIN OBSERVATION t → v(t) RELATED TO THE LINEAR (UPPER GRAPH) AND NONLINEAR (LOWER GRAPH) CASES FOR r = 100%

FIGURE 9 . 5 )Figure 9

 959 FIGURE 9. SENSITIVITY ANALYSIS WITH RESPECT TO THE LOAD RATE;

FIGURE 10 .

 10 FIGURE 10. PROBABILITY DISTRIBUTION FUNCTION OF V ∞ FOR δ = 0.02 (SOLID LINE) AND δ = 0.1 (DASHED LINE)

  display

FIGURE 11 . 2 FIGURE 12 .

 11212 FIGURE 11. CONFIDENCE REGION OF THE DYNAMIC MAG-NIFICATION FACTOR ν → Y LIN (2πν) FOR THE LINEAR CASE WITH δ = 0.02.

FIGURE 13 .

 13 FIGURE 13. CONFIDENCE REGION OF THE DYNAMIC MAG-NIFICATION FACTOR ν → Y NONLIN (2πν) FOR THE GEOMETRIC NONLINEAR CASE WITH δ = 0.02

FIGURE 14 .

 14 FIGURE 14. CONFIDENCE REGION OF THE DYNAMIC MAG-NIFICATION FACTOR ν → Y NONLIN (2πν) FOR THE GEOMETRIC NONLINEAR CASE WITH δ = 0.1

TABLE 1 .

 1 DATA FOR THE FE MODEL

		FE	Nodes DOF
	Sector	12	46	138
	Bladed disk 288 864	2592


Liang Ma
email: liangma@tsinghua.edu.cn

Wei Zhang

Bo Hu

Damien Chablat

Fouad BENNIS b,c , Franc ¸ois Guillaume

Determination of subject-specific muscle fatigue rates under static fatiguing operations

Keywords: static muscular strength, joint strength decline, muscle fatigue rate, subject-specific fatigue rate

Cumulative local muscle fatigue may lead to potential musculoskeletal disorder (MSD) risks , and subject-specific muscle fatigability needs to be considered to reduce potential MSD risks. This study was conducted to determine local muscle fatigue rate at shoulder joint level based on an exponential function derived from a muscle fatigue model. Forty male subjects participated in a fatiguing operation under a static posture with a range of relative force levels (14% -33%). Remaining maximum muscle strengths were measured after different fatiguing sessions. The time course of strength decline was fitted to the exponential function. Subjectspecific fatigue rates of shoulder joint moment strength were determined. Good correspondence (R 2 > 0.8) was found in the regression of the majority (35 out of 40 subjects). Substantial inter-individual variability in fatigue rate was found and discussed.

Introduction

Human intervention is often involved in occupational activities, especially in material handling, assembly, and maintenance tasks (Melhorn et al., 2001a,b;[START_REF] Kumar | Theories of musculoskeletal injury causation[END_REF]. In those activities, the operator needs sufficient muscle strength to meet force requirement for operating equipment or sustaining external loads. Insufficient strength can lead to overexertion of the musculoskeletal system and to consequent injuries [START_REF] Armstrong | A conceptual model for work-related neck and upper-limb musculoskeletal disorders[END_REF][START_REF] Chaffin | Occupational biomechanics, 3rd Edition[END_REF]. R1Q1 A decrease in muscle strength is often experienced in a physical operation under a sub-maximal force, either in a continuous way or in an intermittent way [START_REF] Wood | Minimizing fatigue during repetitive jobs: optimal work-rest schedules[END_REF][START_REF] Yassierli | Effects of age, gender, and task parameters on fatigue development during intermittent isokinetic torso extensions[END_REF]. This decrease in maximum force output results from different sources, such as muscle fatigue, musculoskeletal disorders, lack of motivation, etc. Among those sources, muscle fatigue is one of the most prevalent reasons and is defined as "any exercise-induced reduction in the capacity to generate force or power output" [START_REF] Vøllestad | Measurement of human muscle fatigue[END_REF]. Muscle fatigue exposes operator to more risks of overexertion, and cumulative muscle fatigue may result in musculoskeletal disorders (MSDs) [START_REF] Armstrong | A conceptual model for work-related neck and upper-limb musculoskeletal disorders[END_REF][START_REF] Chaffin | Occupational biomechanics, 3rd Edition[END_REF]. R1Q2 Fatigue progression is closely dependent on task assignment and subject-specific fatigue attributes. Different task parameters (load, duration of exertion, posture and motion, etc.) lead to different fatigue progressions in physical operations [START_REF] Chaffin | The evolving role of biomechanics in prevention of overexertion injuries[END_REF][START_REF] Yassierli | Effects of age, gender, and task parameters on fatigue development during intermittent isokinetic torso extensions[END_REF][START_REF] Enoka | Muscle fatigue-from motor units to clinical symptoms[END_REF]. Individual physical attributes (e.g., strength, fatigue rates, recovery rates, etc.) can influence the fatigue progression as well. It is believed that fatigue attributes differentiate from each other among operators [START_REF] Yassierli | The influence of age on isometric endurance and fatigue is muscle dependent: a study of shoulder abduction and torso extension[END_REF][START_REF] Yassierli | Effects of age, gender, and task parameters on fatigue development during intermittent isokinetic torso extensions[END_REF][START_REF] Avin | Sex differences in fatigue resistance are muscle group dependent[END_REF][START_REF] Avin | Age-related differences in muscle fatigue vary by contraction type: a meta-analysis[END_REF]. Determination of subject-specific fatigue attributes is of interest and of importance for physical work design [START_REF] Chaffin | The evolving role of biomechanics in prevention of overexertion injuries[END_REF].

Muscle fatigue progression has been studied mainly from two different approaches. One approach is maximum endurance time (MET) approach. MET can assess the ability to resist fatigue by measuring the maximal duration while exerting a force at a specific level until failure. A large amount of effort has been contributed to developing MET models for different muscle groups under different static working conditions [START_REF] Rohmert | Ermittlung von Erholungspausen für statische Arbeit des Menschen[END_REF][START_REF] Rohmert | Problems in determining rest allowances Part 1: use of modern methods to evaluate stress and strain in static muscular work[END_REF][START_REF] Rohmert | A study stressing the need for a static postural force model for work analysis[END_REF][START_REF] Bishu | Force-endurance relationship: does it matter if gloves are donned?[END_REF][START_REF] Kanemura | Effect of visual feedback on muscle endurance in normal subjects[END_REF][START_REF] Mathiassen | Prediction of shoulder flexion endurance from personal factor[END_REF][START_REF] Garg | The effect of maximum voluntary contraction on endurance times for the shoulder girdle[END_REF][START_REF] Law | Endurance time is joint-specific: A modelling and metaanalysis investigation[END_REF]. Although the MET models can predict the maximum endurance time under a given relative force level, the decrease in the muscular strength cannot be predicted directly by MET models. Moreover, due to the nature of the formation in those MET models from group data, it is difficult to determine subjectspecific fatigue attributes.

Another approach to characterise muscle fatigue progression is to develop muscle fatigue models, and hence to predict the strength decline directly. Some work [START_REF] Giat | A musculotendon model of the fatigue profiles of paralyzed quadriceps muscle under fes[END_REF][START_REF] Ding | A predictive model of fatigue in human skeletal muscles[END_REF] contributed to complex physiological muscle fatigue models. Those models are able to describe the muscle fatigue progression precisely for a single muscle. However, they are too complex for industrial application due to the difficulty of identifying a great number of parameters in the model. Some researchers established some fatigue models by conducting fatiguing tasks [START_REF] Deeb | An exponential model of isometric muscular fatigue as a function of age and muscle groups[END_REF][START_REF] Søgaard | The effect of sustained low-intensity contractions on supraspinal fatigue in human elbow flexor muscles[END_REF][START_REF] Roman-Liu | Quantitative assessment of upper limb muscle fatigue depending on the conditions of repetitive task load[END_REF], 2005;[START_REF] Iridiastadi | Muscle fatigue and endurance during repetitive intermittent static efforts: development of prediction models[END_REF]. Among those models, exponential declines in muscular strength are found [START_REF] Deeb | An exponential model of isometric muscular fatigue as a function of age and muscle groups[END_REF][START_REF] Roman-Liu | Quantitative assessment of upper limb muscle fatigue depending on the conditions of repetitive task load[END_REF][START_REF] Yassierli | The influence of age on isometric endurance and fatigue is muscle dependent: a study of shoulder abduction and torso extension[END_REF]. However, the fitting param-R2Q9 eters in those exponential functions could not implicate more information about fatigue attributes of each subject. Some other researchers [START_REF] Liu | A dynamical model of muscle activation, fatigue, and recovery[END_REF][START_REF] Ma | A new simple dynamic muscle fatigue model and its validation[END_REF][START_REF] Xia | A theoretical approach for modeling peripheral muscle fatigue and recovery[END_REF] have tried other models to describe muscle fatigue progression. Xia and Frey Law R1Q3 (2008) proposed a three-compartment muscle fatigue model based on muscle motor units model in [START_REF] Liu | A dynamical model of muscle activation, fatigue, and recovery[END_REF], and they have run simulation to demonstrate fatigue progression under a variety of loading conditions. In this model, fatigue and recovery attributes of different types of muscle fibers were assigned with different values in the simulation. However, the lack of validation limits the application of this model. [START_REF] Ma | A new simple dynamic muscle fatigue model and its validation[END_REF] proposed a muscle fatigue model to describe muscle fatigue progression from a macro perspective. In this model, task parameters and muscle fatigue rate are combined together to understand the fatigue caused by tasks and fatigue attributes. [START_REF] Ma | A novel approach for determining fatigue resistances of different muscle groups in static cases[END_REF] developed an approach to determine fatigue resistances of different muscle groups using this fatigue model. Twenty-four MET models [START_REF] El Ahrache | Percentile values for determining maximum endurance times for static muscular work[END_REF] for different muscle groups can be effectively fitted and explained by this approach. We suggest that this muscle fatigue model is ca-R2Q10 pable of assessing fatigue progression of a muscle group in static cases. Moreover, we found that the muscle fatigue progression of each subject under static cases can be predicted in the form of an exponential function derived from the fatigue model, and one important factor (fatigue rate) in this model emerges to represent subject-specific muscle fatigue attribute.

Regarding the subject-specific fatigability, some other measures were used to R1Q0 assess muscle fatigability, such as endurance time, EMG power spectrum (Median frequency (MF) and Median Power Frequency (MPF)), Mean Arterial Pressure (MAP) and so on [START_REF] Clark | Gender differences in skeletal muscle fatigability are related to contraction type and emg spectral compression[END_REF][START_REF] Hunter | Fatigability of the elbow flexor muscles for a sustained submaximal contraction is similar in men and women matched for strength[END_REF][START_REF] Hunter | Muscle endurance is greater for old men compared with strength-matched young men[END_REF][START_REF] Yoon | Mechanisms of fatigue differ after low-and high-force fatiguing contractions in men and women[END_REF][START_REF] Frey Law | Endurance time is joint-specific: a modelling and meta-analysis investigation[END_REF][START_REF] Côté | A critical review on physical factors and functional characteristics that may explain a sex/gender difference in work-related neck/shoulder disorders[END_REF]. However, as pointed out by [START_REF] Vøllestad | Measurement of human muscle fatigue[END_REF], the greatest limitation is that those measures are indirect to measure muscle fatigue. Therefore, we chose to use the fatigue rate in [START_REF] Ma | A novel approach for determining fatigue resistances of different muscle groups in static cases[END_REF] to represent inter-individual difference in fatigue progression to beyond those limitations.

We conducted this study to verify whether the fatigue progression under a static R1Q4 operation can be well fitted by a specific exponential function derived from the R2Q11 muscle fatigue model and to check whether the fitting parameter (fatigue rate) could represent subject-specific fatigue attributes among different subjects. This paper is organised as follows: Section 2 describes the theoretical approach to determine subject-specific fatigue rate based on a muscle fatigue model. Section 3 presents materials, methods, and experiment procedure in this study. Section 4 and Section 5 show results and provide discussion.

R2Q12, R1Q17

Subject-specific fatigue rate determination

The purpose of this study was to determine subject-specific muscle group fatigue attributes by analysing muscle fatigue progression based on a theoretical muscle fatigue model. [START_REF] Ma | A new simple dynamic muscle fatigue model and its validation[END_REF] proposed a muscle fatigue model in the form of a differential equation (Eq. 1). The muscle fatigue model describes the change of the maximum R1Q6,

R2Q13 R2Q14

remaining strength over time. Related parameters and their descriptions are given in Table 1. In this model, the fatigue rate (k) is a parameter to indicate the relative speed of strength decline within a muscle.

dF rem (t) dt = -k F rem (t) MVC F load (t) (1)
[Table 1 about here.]

In a static muscular operation, F load (t) keeps constant, and the reduction of the R1Q7 muscular strength can be further predicted by Eq. 2. This equation describes the muscle fatigue progression in the form of an exponential function. Three parame-R2Q13 ters (F max , F load , and k) act upon the fatigue progression under a static operation.

In general, F load is determined by the task design, and it can be measured or calculated via force analysis, and F max and F rem can be measured to unfold the muscle fatigue progression. The fatigue rate k is task independent and it is influenced by several factors (e.g., muscle fiber type composition, age, and gender) [START_REF] Ma | A novel approach for determining fatigue resistances of different muscle groups in static cases[END_REF]. R1Q12

F rem (t) F max = e -k f MVC t
(2) R1Q8 According to the definition of muscle fatigue, muscle fatigue progression can also be described by measuring the maximum remaining muscle strengths over time during a fatiguing operation. If the same muscle progression can be depicted using both ways, it will suggest the parameter k could be determined using the muscle fatigue model. Therefore, it is essential to verify whether the fatigue progression of each subject under a static fatiguing operation follows an exponential function in the form of Eq. 2 with a high coefficient of determination R 2 . R1Q9, R1Q10,

R1Q11,R2Q6

Suppose that we have already a set of real measurements for a given static operation, where F t i indicates the maximum remaining strength F rem at time instant t i . At the beginning of a physical task, the subject is supposed having no muscle fatigue. Therefore, the F t=0 can be treated as the maximum voluntary contraction F max . Equation 2 can be further transformed into Eq. 3.

ln

F t i MVC f MVC = -k t i (3)
Since f load and MVC are measured and/or known, fatigue rate k can be further determined by linear regression. A high goodness of fit between maximum remaining strengths and the exponential function would suggest the usefulness of the fatigue rate. The goodness of fit is assessed by the R 2 value in a linear regression without an intercept.

Materials and Methods

Subjects

Since the focus of this study is on manufacturing and assembly and the majority of the operators are male workers, 40 right-handed male industrial workers participated in the experiment after signing a written informed consent. The age, stature, body mass, upper limb anthropometry data, and body mass index (BMI)

were recorded or measured upon arrival at the laboratory (see Table 2). Participation was limited to individuals with no previous history of upper limb problems.

Ethical approval for this study was obtained from the human research ethical advisory committee of Tsinghua University.

[Table 2 about here.]

Task design

In this study, we used a typical overhead drilling operation under laboratory conditions to measure muscle fatigue progression at shoulder joint level (see Fig. 1). This task was simplified from a real drilling operation from the program of the European Aeronautic Defence and Space (EADS). This operation was selected as a typical task because there are a few ergonomics problems (Melhorn et al., 2001a). R2Q15 A heavy external load demands great physical strength to hold a machine and maintain the operation for a certain period, and local muscle fatigue occurs rapidly in upper limb and lower back. MSD risks can be increased by force overexertion and long-lasting vibration while drilling.

The magnitude of the external load and the duration of the operation are two key factors to simulate the real situation. The relative load needs to be adequate so that subjects can experience fatigue and can endure the operation for a certain period. In this case, according to the strength model in [START_REF] Chaffin | Occupational biomechanics, 3rd Edition[END_REF] and the MET models in El ahrache et al. ( 2006), a subject must apply a drilling force of 25 N and hold a drilling machine with a mass of 2.5 kg. The drilling force is only applied along the drilling direction towards the subject. The estimated moment generated by the external load (including the weight of the arm) is about 33% of the shoulder joint flexion strength of a 50 th percentile male and the endurance time under this load is estimated around 4 minutes.

[Figure 1 about here.]

Measures

In this study, we used shoulder joint moment strength to describe fatigue progression and measured the maximum force output in the drilling direction to estimate shoulder joint moment strength (see Fig. 2). R1Q5 We assumed that the measured force was determined by joint moment strength of the right upper limb. Shoulder joint and elbow joint have similar strength profiles according to the joint moment strength models [START_REF] Chaffin | Occupational biomechanics, 3rd Edition[END_REF], and shoulder joint has higher fatigability in MET models than elbow joint (Frey [START_REF] Law | Endurance time is joint-specific: A modelling and metaanalysis investigation[END_REF]. Furthermore, it was obvious that shoulder joint was charged with a much larger moment load than elbow joint in this drilling case. Therefore, the bottleneck for the output strength was shoulder joint.

The moment load in shoulder joint can be approximately estimated by Eq. 4 (see Fig. 2).The mass and the centre of gravity of each body segment were estimated from the anthropometry database [START_REF] Chaffin | Occupational biomechanics, 3rd Edition[END_REF].

Γ load = s -e 2 × G u + w + e 2 -s × G f + d + w 2 -s × G m + (d -s) × F d (4)
where s, e, w, and d represent the coordinates of the positioning sensors attached to the shoulder (S), elbow (E), wrist (W), and drilling contact point (D), respectively.

In our experiment, since the subject's arm is strictly limited within the sagittal plane, we just measured the force and calculated the torque within this plane. R1Q16

[Figure 2 about here.]

Since we use joint moment strength, the fatigue model in Eq. 1 can be changed to Eq. 5 by replacing all of the force terms with joint moment terms.

dΓ rem (t) dt = -k Γ load (t) Γ max Γ rem (t) (5)
Under this simplified case, Γ load can be determined by force analysis, and f MVC can be calculated from data analysis by normalizing the Γ load over Γ max . Γ t i is a joint maximum remaining strength, and it can be estimated by measuring the maximum remaining force F t i in the drilling direction.

Material

In the experiment, we want to measure the muscular strength in the drilling direction and to estimate the joint moment strength around the shoulder joint. Therefore, force measurement and motion capture devices were used (see Figure 1).

We used a dynamometer to measure the drilling force in the drilling direction (see Figure 1). The dynamometer measures the pressing force perpendicular to the load cell surface with a measurement range upto 300 N and a precision of 1 N. R1Q13 We use the magnetic motion capture device FASTRAK R (POLHEMUS Inc.) to capture the upper limb posture in the experiment. As shown in Figure 1, we attached four positioning sensors to the key joints of upper limb and the drilling machine. We captured the Cartesian coordinates of the shoulder, the elbow, the wrist, and the contact point between the drilling machine and the work piece. The tracking device runs at 30 Hz per sensor and has a static position accuracy of 1 mm. We used the recorded coordinates of each tracker to reconstruct the posture of the worker in post-experiment analysis.

We provided the drilling force with a wooden beam with a mass of 10 kg.

We used wooden material to avoid magnetic distortions caused by ferrous material and to ensure motion capture accuracy. We suspended the wooden beam with an inclination angle between the beam and the horizontal line of 14.5 • . During operation, the subject had to push the beam against the force measurement device and hold it for a certain period. According to the force analysis of the pendulum, a tangential force of 25 N was charged to the upper limb. Before each trial, we calibrated this external load to ensure that there was exactly a force of 25 N applied in the drilling direction. We provided the weight of the drilling machine with a drilling tool made from concrete with a mass of 2.5 kg.

Experiment Procedure

Each subject had to complete ten sessions: one MVC session and nine fatiguing sessions.

In the MVC session, maximum voluntary contraction (MVC) was determined as the greatest exerted force in the drilling direction over three trials. In each trial, we verbally encouraged subjects to maintain the maximum force peak for three to five seconds. The measured MVC was also denoted as F 0 to represent the subject's initial maximum strength at the beginning of the operation. Between each trial, subjects were asked to take at least a 5-minute rest until self-reported full recovery [START_REF] Chaffin | Occupational biomechanics, 3rd Edition[END_REF].

There were nine fatiguing sessions with different time intervals (15, 30, 45, 60, 75, 90, 120, 150, and 180 seconds). The sequence to complete those nine sessions was randomly assigned for each subject. In each fatiguing session, subject was asked to hold the constant external load for the time interval t i (e.g., 30 sec). After that, the remaining muscle strength F t i (e.g., F 30 ) was measured by asking the subject to exert the maximal voluntary strength with a force peak from three to five seconds. After the measurement, subjects took a rest for at least five minutes or even longer until self-reported total recovery. After the recovery, the subject was asked randomly to conduct a MVC trial. R2Q2 Full recovery would be recognized if the measured maximum strength in this trial was more than 95% of the MVC. Otherwise, the subject would be asked to take longer break until full recovery. Once subject reported that he could not sustain the R1Q15 operation within the session, experiment would be stopped immediately to avoid injuries to subject.

Within each session, subject was seated upright, and right shoulder was fixed to a shoulder support against the wall to restrict the movement of shoulder and to decrease the engagement of lower back. Left upper limb was free, and right upper limb was limited in the sagittal plane by position constraints. The position constraints provided posture references to subject to maintain the initial posture as well as possible, but provided no support to upper limb. The seated height and location was adjusted according to subjects height and upper limb length to reduce variability among the different subjects. R1Q14

Data Analysis

The objective of this analysis was twofold. R2Q16

1. to test if the fatigue progression of each subject in shoulder joint maximum strength can be well fitted by the exponential function derived from the muscle fatigue model or not; R2Q3 2. to analyse the relationship between muscle fatigue rate and joint moment strength.

For the first objective, we fitted muscle fatigue progression of each subject with the exponential function (Eq. 3). The coefficient of determination of each fitting was recorded for each subject and analysed. For the second objective, we selected two groups of subjects to assess the relationship of joint moment strength and joint fatigue rate, since muscles engaged in the action mainly determine maximum joint strength and we assumed that determinant muscle-related factors for muscle strength could probably act effects on muscle fatigue rate as well. One group (Group A) is the subjects with 10 highest joint moment strengths; another group is the subjects with 10 lowest joint moment strengths. Besides joint strength, some other measures (e.g., BMI, age) may also influence fatigue rate. Correlations between fatigue rate and other measures were also statistically analysed. We used SPSS to do all the statistical analysis and fitting. R1Q12, R1Q17, R1Q19 4. Results

Fatigue progression

Joint moments of each subject were normalized over the estimated maximum moment strength. Then the normalized values were fitted using Eq. 3 to calculate the fitting coefficient R 2 and to determine the fatigue rate. The statistical results of R 2 of the regression and the fatigue rate k were listed in Table 3, and the histograms of R 2 and k were shown in Fig. 4 and Fig. 5.

35 out of 40 subjects had a high coefficient of determination R 2 over 0.8 in joint moment regression. Four out of the other five subjects had a fair coefficient R 2 over 0.63, and only one of the five subjects had a very poor R 2 (0.23). R2Q4

[Table 3 about here.]

[Figure 3 about here.]

Discussion

Muscle fatigue progression

In this study, we found that the fatigue progression at shoulder joint among most of the subjects (35/40) can be well fitted (R 2 > 0.8) by the muscle fatigue model. Five out of 40 subjects were found with low R 2 coefficients. R2Q6 There are some reasons leading to those relatively poor fittings. First, the motivation and the attitude of the subject during the experiment could influence the result. Second, even though the posture of the arm was strictly constrained in the R2Q17 sagittal plane, subjects could still have a certain degree of mobility. The willingness to maintain the posture may influence the muscle recruitment strategy, muscle coordination, and the posture during the experiment. Third, weak muscular strength could probably lead to poor fitting performance indirectly. The subject R2Q20 with a R 2 = 0.23 had the second lowest MVC among the 40 subjects. Lower muscle strength means relatively higher physical workload during the experiment and higher demand to maintain the posture, and hence the static operation would be more possibly aversely influenced due to the unwillingness to maintain the posture.

R2Q7,

R2Q19

Fatigue may occur at any step along the pathway that is involved in muscle contraction [START_REF] Berne | Physiology[END_REF]. Both the metabolic factors within the muscle and the impairment of activation could contribute to the decline in muscle power output [START_REF] Chaffin | Occupational biomechanics, 3rd Edition[END_REF][START_REF] Allen | Skeletal muscle fatigue: cellular mechanisms[END_REF]. In a physical operation, muscle fatigue progression could be influenced by physical task, motivation, and individual fatigue attributes. However, under static operation and effective verbal encour-R1Q21 agement, the influences from motivation could be rather limited. Therefore, muscle fatigue progression was probably mainly caused by relative loads and subjectspecific fatigue attributes in this study.

The fatigue model (Eq. 2) is formed from a macro perspective and can be explained based on motor units principle from a micro perspective [START_REF] Ma | A new simple dynamic muscle fatigue model and its validation[END_REF].The product of relative load and fatigue rate determines the decline of muscle strength.

According to muscle physiology, fatigue rate can be recognized as a parameter to R1Q23 represent the overall fatigue resistant performance of a muscle group at joint level under a specific task [START_REF] Ma | A novel approach for determining fatigue resistances of different muscle groups in static cases[END_REF].

We selected this model to determine fatigue attribute for the following reasons: (1) in comparison to Deeb's model [START_REF] Deeb | An exponential model of isometric muscular fatigue as a function of age and muscle groups[END_REF], this model enables us to decouple relative load and subject-specific attribute; (2) muscle fatigue rate can be influenced by muscle composition, neuromuscular activation patterns, and coordination among single muscles, therefore the fatigue rate in this model could cover more effects of influencing factors than the fatigue rates of different types of muscle fibers in [START_REF] Xia | A theoretical approach for modeling peripheral muscle fatigue and recovery[END_REF]; (3)this model is less complex than the R1Q22 three-compartment model in [START_REF] Xia | A theoretical approach for modeling peripheral muscle fatigue and recovery[END_REF] , and it could be practical for industry application.

Subject-specific fatigue rate

We found that there are substantial differences among different subjects in fatigue rates at shoulder joint level. It suggests that the fatigue rate k can be used to characterise different fatigue attributes among subjects at shoulder joint level under this specific fatiguing condition.

The underlying mechanisms for those differences in fatigue rate are very complex. Physiologically, differences in fatigue rates are mainly caused by four factors: (1) muscle strength (muscle mass) and associated vascular occlusion, (2) substrate utilisation, (3) muscle composition, and (4) neuromuscular activation patterns [START_REF] Hicks | Sex differences in human skeletal muscle fatigue[END_REF][START_REF] Ma | A novel approach for determining fatigue resistances of different muscle groups in static cases[END_REF]. At the same time, demographic parameters (age, gender) and their interactions can lead to changes in muscle composition as well [START_REF] Mademli | Effect of voluntary activation on age-related muscle fatigue resistance[END_REF][START_REF] Yassierli | Effects of age, gender, and task parameters on fatigue development during intermittent isokinetic torso extensions[END_REF]. In addition, personal working experience or physical exercises and living style (e.g., smoking) can also change the muscle strength and endurance via the adaptive response of muscle cells to regular external loads [START_REF] Berne | Physiology[END_REF][START_REF] Wüst | Skeletal muscle properties and fatigue resistance in relation to smoking history[END_REF]. All those factors generate effects together on subject-specific fatigue attributes.

R1Q23 We found also that fatigue rates are positively correlated to maximum joint moment strength in this study, even though the relative loads for the subjects with higher strength were lower. Between-subject differences in the ratio of type I mus-R2Q17, R2Q21 cle fibers (slow twitch, more fatigue resistant) to type II muscle fibers (fast twitch, less fatigue resistant) might explain the differences in fatigue rates. Muscle strength depends strongly on muscle fiber size and muscle fiber composition [START_REF] Fitts | The determinants of skeletal muscle force and power: their adaptability with changes in activity pattern[END_REF]. Subjects in Group A and Group B did not have significant differences in BMI and age, which implicates that the strength differences were probably not mainly caused by muscle mass or muscle fiber size or age-related factors. It could be concluded that the strength differences were caused by different compositions of muscle fiber types. Subjects with higher strengths could probably have a higher proportion of type II muscle fibers and lower proportion of type I muscle fibers.

That leads to higher joint moment strength and faster fatigability in the muscle.

We did not find significant correlation between fatigue rate and age and between fatigue rate and BMI. Regarding the age effect, we did not control our subjects strictly to two age groups. The subjects in this experiment were in their young age or middle age, which might not enough to reveal the aging effect. Regarding the BMI, most of the subjects were in normal weight and overweight group. Only a few subjects belong to Class I obesity or underweight category.

Posture change

Although posture references were provided to avoid mismatches in different test trials, it was still very difficult for subjects to maintain the purely static posture during the operation. Small changes occurred in the experiment, but those changes did not generate excessive variation in joint strength. In our case, the variation of the maximum joint moment strength is no more than 3% relative to the maximum strength under the initial posture according to the joint moment strength model [START_REF] Chaffin | Occupational biomechanics, 3rd Edition[END_REF]. The change of the joint strength due to posture change might lead change of the relative strength. We checked the sensitivity of the change, and the changed maximum joint strength would lead to no more than4% change of the relative strength (4% f MVC ), which was acceptable in this case. R1Q20 Several reasons may cause posture change in the experiment. Fatigue might be one of those reasons. Changes in the posture can be explained by a global posture control strategy, which includes decreasing the joint loads in the operation by moving the upper limb closer to the body; a similar finding has been reported by [START_REF] Fuller | Posture-movement changes following repetitive motion-induced shoulder muscle fatigue[END_REF]. However, that change would influence joint strength(Roman- Liu and Tokarski, 2005;[START_REF] Anderson | Maximum voluntary joint torque as a function of joint angle and angular velocity: Model development and application to the lower limb[END_REF]

R2Q23

There are several limitations in this study. First, the focus of the present study is on the fatigue effect in static industrial operations in a continuous working process.

Recovery effect is not taken into consideration in this study. Second, this study was conducted under a simplified overhead drilling operation, and the conclusion drawn from this study has a strong task dependency. A simplified overhead drilling operation decreases the reliability of applying the findings into industry. Some other MSD causes, such as vibration [START_REF] Kattel | The effects of rivet guns on hand-arm vibration[END_REF] were neglected from this study. Third, the force analysis in this study is available only for a static case.

In a real operation, the motion involved in the operation could result in a different dynamic workload. Moreover, only fatigue with the relative force falling from 14% to 33% (Mean=24.3%, SD=4.4%) of the specific job operation was tested, so the result that was obtained is available only for similar physical operations. Last but not least, the fatigue progression was measured under static isometric contraction, and the results could not be extended for dynamic or intermittent fatiguing tasks.

Conclusions and perspectives

This paper provides an experimental approach to determine subject-specific fatigue rate at shoulder joint level. Fatigue progression in a simplified static drilling operation was measured and analysed using an exponential muscle fatigue model.

q 2 S E W D G u G f G m F d

. Relationship between fatigue rate and joint maximum strength

A pair-wise correlation matrix was determined between joint moment fatigue rate, joint moment strength, age and BMI. The results were shown in Table 4. It was found that joint moment strength is strongly correlated with fatigue rate (p < 0.05), while no strong correlations were found in the pair of fatigue rate and BMI and the pair of fatigue rate and age.

R1Q19

No significant differences were found between two groups in age (Group A: 42.9 (SD=7.4); Group B: 39.1 (SD=15.4), p = 0.49) and BMI (Group A: 24.9 (SD=2.4); Group B: 25.2 (SD=4.2) p = 0.78). Significant differences (p = 0.00) were found in joint maximum strength between Group A (mean=60.8 Nm (SD=4.7Nm)) and Group B (Mean=37.7 Nm (SD=3.3 Nm)) . R1Q20 With use of t-test, Table 5 showed the difference of fatigue rates between different groups. The subjects with higher strength have significantly higher fatigue rate even though the relative load is smaller than the subjects with lower strength. R2Q21 We calculated the posture of upper limb during the drilling operation from the motion data. Because the arm was constrained in the sagittal plane, only the flexion angles of shoulder joint and elbow joint were calculated to represent arm posture to eliminate influence from different limb lengths. The statistical results of elbow

  . Besides fatigue, there were still other error sources leading to the change of the posture. First, the actual posture was determined by the anthropometry of different subjects. Different arm lengths could cause potential differences in elbow flexion and shoulder flexion. Second, the posture was calculated from the position sensors attached to the key joints. Each subject might have different sensor configurations, which might lead to calculation errors. Third, there might be differences among the postures that each subject took in

Figure 1 :

 1 Figure 1: Seated static posture in the experiment and materials used in the experiment. F drill : drilling force; G machine : the weight of the drilling machine; G upper : the weight of the upper arm; G lower : the weight of the lower arm

Figure 2 :Figure 3 :

 23 Figure 2: Force schema in the drilling operation. S: shoulder; E: elbow; W: wrist; D: drilling point. F d : drilling force applied to the upper limb; G m : weight of the drilling machine; G u : weight of the upper arm; G f : weight of the forearm; q 1 : flexion of the upper arm; and q 2 : flexion of the elbow

Figure 4 :

 4 Figure 4: Histogram of coefficient of determination R 2 R1Q18 R1Q18 6.1. Relationship between fatigue rate and joint maximum strength

Figure 5 :

 5 Figure 5: Histogram of fatigue rate at shoulder joint

Table 1 :

 1 Parameters in the dynamic fatigue model

	Item	Unit Description
	MVC or F max	N	Maximum voluntary muscle strength under nonfatigued
			state
	F rem (t)		Maximum voluntary remaining muscle strength at time t
	F load (t)	N	External load that the muscle needs to bear
	k	min -1 Fatigue rate
	%MVC		Percentage of the voluntary maximum contraction
	f MVC		%MVC/100, f MVC = F load MVC .

Table 3 :

 3 Statistical analysis of fatigue rate k

			Item Mean SD Minimum Maximum
	Joint moment strength (Nm) 45.1	7.4	67.4	32.1
				k 1.02 0.49	0.37	2.29
			R 2 0.87 0.14	0.23	0.99
						R2Q18
		20
	Number of subjects	10
		0
		0.2	0.4	0.6	0.8	1.0
				R 2

Table 4 :

 4 Correlation matrix for study variables.( * p < 0.05)

		Fatigue rate BMI Joint moment strength Age
	Fatigue rate	1	-0.09 0.616 *	0.033
	BMI		1	0.09	0.40 *
	Joint moment strength			1	0.072
	Age				1
			32

Table 5 :

 5 Effect of muscle strength on fatigue rate

	Group Mean SD	t	p-Value
	k
	A	1.47	0.53 4.628 0.0001
	B	0.64	0.20
				R2Q18
	6.2. Posture change during the drilling operation

Acknowledgments

This research was supported by the EADS and the Région des Pays de la Loire (France) in the context of collaboration between the École Centrale de Nantes (Nantes, France) and Tsinghua University (Beijing, PR China). This research was supported by the National Natural Science Foundation of China as well(Grant number: NSFC71101079). R2Q24

flexion and shoulder flexion angles across participants were calculated and shown in Table 6. The posture change during the working process is shown in Figure 3. The changes in the posture followed the following tendency: the greater the fatigue was, the closer the upper limb was to the trunk. In this way, the moment produced by the mass of the upper limb around shoulder joint could be reduced.