A Batou
email: anas.batou@liverpool.ac.uk

A sensitivity-based one-parameter-at-a-time model updating method

Keywords: Model Updating, Sobol Indices, Sensitivity Indices, Inverse problems )

This paper is interested in model updating problems which consists in identifying optimal values of model parameters by comparing the model outputs with the experimental outputs. Such a problem generally yields a challenging multivariate inverse problem to be solved in high dimension. The highdimensionality requires the use of a global optimization algorithm in order the explore efficiently the parameters space. In this paper we propose an alternative algorithm which allows each model parameters to be identified separately and sequentially by solving separated univariate inverse problems. For each parameter, a devoted inverse problem is introduced by identifying an output which is sensitive to this parameter only, the sensitivity being quantified using Sobol indices. The proposed method is illustrated through a three-storey structure for which experimental measurements are collected.

Introduction

In the context of model validation, the model updating methods consist in calibrating a computational model in order to reduce the distance between the experimental outputs and the outputs provided by the computational model [START_REF] Friswell | Finite Element Model Updating in Structural Dynamics[END_REF]. The model updating methods of computational models have been intensively studied during the last four decades. In structural dynamics, efficient methods have been proposed (see for instance [START_REF] Ljung | System Identification: Theory for the User[END_REF][START_REF] Walter | Identification of Parametric Models from Experimental Data[END_REF][START_REF] Friswell | Finite Element Model Updating in Structural Dynamics[END_REF]) and are now commonly used in industry. There are two main types of methods: the direct methods (see for instance [START_REF] Baruch | Optimisation procedure to correct stiffness and flexibility matrices using vibration data[END_REF][START_REF] Friswell | Direct updating of damping and stiffness matrices[END_REF][START_REF] Schultz | Model updating using constrained eigenstructure assignment[END_REF][START_REF] Zimmerman | Model correlation using multiple static load and vibration tests[END_REF]) which consist in directly modifying the stiffness and mass matrices and the indirect methods which consists in updating some physical parameters (see for instance [START_REF] Bollinger | Finite element model updating for FEA/EMA modal correlation via constrained optimization theory[END_REF][START_REF] Cogan | An improved frequency response residual for model correction[END_REF][START_REF] Jaishi | Structural Finite Element Model Updating Using Ambient Vibration Test Results[END_REF]) of the model. The latter method can be described in three steps: [START_REF] Baruch | Optimisation procedure to correct stiffness and flexibility matrices using vibration data[END_REF] The first one consists in constructing a computational model for which the parameters are set to nominal values. These parameters can be related to material properties (Young's modulus, mass density, and so on), geometry (CAD, thickness, area moments of inertia, and so on) and boundary conditions. [START_REF] Baruch | Modal data are insufficient for identification of both mass and stiffness matrices[END_REF] The second step consists in selecting the parameters to be updated. This step generally relies on engineering expertise combined with a sensitivity analysis of the quantities of interest with respect to the candidate parameters to be updated (see for instance [START_REF] Ben-Haim | Selective sensitivity in the frequency domain: Theory[END_REF][START_REF] Cogan | Model order reduction by selective sensitivity[END_REF][START_REF] Mottershead | The sensitivity method in finite element model updating: A tutorial[END_REF]). In [START_REF] Yuan | Parameter selection for model updating with global sensitivity analysis[END_REF], for updating stochastic models, the authors have performed the model selection by defining a composite sensitivity index that discriminates explicitly between sets of parameters with correctly-modelled and erroneous statistics. [START_REF] Ben-Haim | Selective sensitivity in the frequency domain: Theory[END_REF] The third step consists in updating the mostsensitive parameters using experimental data. In general, this step is carried out by reducing a "distance" between the experimental outputs and the corresponding outputs computed using the computational model.

The number of most-sensitive parameters is often very large and the updating of those later generally requires the solving of a non-convex inverse problem in high dimension. This high dimension prohibits the use of grid search methods and global optimization methods are required such as genetic algorithm, particle swarm optimization or simulated annealing [START_REF] Pardalos | Handbook of global optimization[END_REF]. These classes of methods require the fitting of solver parameters in order to be efficient and enable the escape from local minima.

The objective of the new approach proposed in the present paper consists in addressing the problem of the high dimension in model updating method by replacing the original inverse problem by a sequence of small dimension (ideally dimension one) inverse problems which can be solved using a grid search approach. To achieve this objective the proposed approach consists in transforming the observed outputs of the model into a set of new outputs. Each new output is sensitive to only one parameter to be updated and is not sensitive to the other parameters. Then by comparing the new outputs of the model to the corresponding experimental ones, it is possible to define a sequence of separated inverse problems in order to update all the parameters one at a time and thus it is possible to use a grid search algorithm solver for each one-dimensional inverse problem.

In Section 2, the new model updating method is presented. In this section, first the construction of the new outputs is presented, then the model updating method using these outputs is presented. The approach is illustrated through simple examples. In Section 3 the updating method is validated experimentally through a 3-storey structure for which experimental data has been collected.

Model updating algorithm

Output transformation

In this paper we consider a computational model for which the output vector y(x) ∈ R m is a function of the vector x ∈ C x ⊂ R n of the model parameters to be updated. It is assumed that a corresponding experimental output y exp has been collected. At this stage, for readability, it is assumed that the fidelity of model is high and therefore, the experimental response can be predicted by the computational model if the true values of the model parameters are used. The presence of modelling errors will be addressed in Section 2.5. A classical way to identify the model parameter vector x consists in finding the value x opt of the vector x which makes the model output y(x opt ) as close as possible to the experimental output y exp . The distance between the two outputs can be measured using the mean square distance y(x opt )y exp for instance, i.e.,

x opt = arg min x∈Cx y(x)y exp .

(

) 1
The solving of the multivariate inverse problem (1) generally requires the use of a global optimization algorithm such as the simulated annealing algorithm, genetic algorithms or particle swarm optimization, in order to explore the parameter space efficiently [START_REF] Pardalos | Handbook of global optimization[END_REF]. The efficiency of these methods depends on the correct tuning of the algorithm parameters which are problem dependent.

A bad tuning of these parameters may yield to a slow convergence of the algorithm and a possible trap in a local minimum. The difficulty of exploring efficiently the parameter space increases with the number of parameter The method proposed in this paper proceeds differently. For each parameter {x i , i = 1, . . . , n}, the output y(x) is transformed into a new output z i (x) with values in R:

z i (x) = g i (y(x)), (2)
which is sensitive to the parameter x i only and thus is not sensitive to the other model parameters. Furthermore it is assumed that this transform is a linear combination of n g basis functions {f j , j = 1, . . . , n g }, i.e,

z i (x) = α i T f (y(x)), (3)
where f = (f 1 , . . . , f ng ) is the vector of the basis functions and α i = (α i,1 , . . . , α i,ng ) is the vector of the coordinates in this basis. In the particular case where f (y(x)) = y(x), the transform g i is just a linear combination of the outputs of the model.

Calculation of the output parameters

For each model parameter x i , the transform g i and thus the vector α i need to be calculated. Since the objective is that each output z i (x) is devoted to the identification of the model parameter x i only, we calculate the vector α i such that the output z i (x) is sensitive to the model parameter x i only. By replacing the model parameter vector x by the random vector X, Sobol indices [START_REF] Sobol | Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates[END_REF][START_REF] Saltelli | Variance based sensitivity analysis of model output. Design and estimator for the total sensitivity index[END_REF] enable the total variance of each random variable X i to be quantified. The advantage of these total-effect sensitivity indices is that they are not local and thus do not require to set values for the other parameters. For each output Z i = z i (X) and for each random variable X j , the total-effect Sobol index Si,j is defined by Si

,j = E X ∼j Var X j (Z i | X ∼j ) Var(Z i ) , (4)
where X ∼j is a random vector for which the components are copies of all the components of X except X j . Then using Eq.( 3) yields

Si,j = α i T [E X ∼j Cov X j (F | X ∼j ) ]α i α i T [Cov(F )]α i , (5)
with

F = f (y(X)). Denoting [C] = [Cov(F )] and [C j ] = [E X ∼j Cov X j (F | X ∼j ) ], the later equation is then rewritten Si,j (α i ) = α i T [C j ]α i α i T [C]α i . (6)
We can see in Eq.( 6) that if matrices [C] and [C j ] are estimated (see [START_REF] Saltelli | Variance based sensitivity analysis of model output. Design and estimator for the total sensitivity index[END_REF] for the estimation of Sobol indices) then the Sobol index Si,j (α i ) can easily be calculated for any value of the vector α i . This will greatly facilitate the calculation of α i .

As explained before, for each parameter x i , the vector α i defining the corresponding output z i (x) has to be calculated such that this output is sensitive to x i only. It can easily be shown that the value α i opt of α i satisfying this condition is such that

Si,i (α i opt ) = 1, Si,j (α i opt ) = 0, for j = i. (7)
Solving Eq.( 7) is equivalent to solving the set of equations

S i,i (α i opt ) = 1, S i,j (α i opt ) = 0, for j = i. (8)
with,

S i,j (α i ) = α i T [C j ]α i . (9)
Indeed if Eq.( 8) is satisfied then

Z i = z i (X) is only sensitive to X i and then (α i opt ) T [C i ]α i opt = (α i opt ) T [C]α i opt , i.e. Si,i (α i opt ) = 1.
The solution of Eq.( 8) can be calculated by solving the following mean square problem,

α i opt = arg min α i ∈R ng Γ(α i ) . (10)
where

Γ(α i ) = (S i,i (α i ) -1) 2 + i =j (S i,j (α i )) 2 + λ α i T α i . (11)
The last term in the right-hand side of Eq.( 11) is a Tikhonov regularization term ( [START_REF] Tikhonov | Solutions of Ill-Posed Problems[END_REF]). The gradient vector ∇Γ(α i ) and the Hessian matrix [H(α i )] of function Γ(α i ) are written as

∇Γ(α i ) = 4 (S i,i (α i ) -1) [C i ]α i + 4 i =j {S i,j (α i ) [C j ]α i } + 2 λ α i , (12)
[

H(α i )] = 4 (S i,i (α i ) -1) [C i ] + 8 [C i ]α i α i T [C i ] T +4 i =j {S i,j (α i ) [C j ] + 2 ([C j ]α i α i T [C j ] T } + 2 λ [I n ]. ( 13
)
Analysing the Hessian matrix [H(α i )], it can be deduced that since S i,i (α i ) ≤ 1, the Hessian matrix is not always positive definite unless λ is sufficiently large. Nevertheless, in the neighbourhood of the solution, S i,i is close to one and therefore, the Hessian matrix is locally positive definite and thus function Γ(α i ) becomes convex when approaching the solution. The regularization parameter λ in Eq.( 11) has to be chosen sufficiently large to guaranty the convexity of Γ(α i ) and sufficiently small to limit the impact of the regularization term. The optimization problem defined by Eqs. ( 10) and ( 11) can be solved using any minimization algorithm. Here we suggest to solve it using the Newton iterative method written as

α i i+1 = α i i -[H(α i i )] -1 ∇Γ(α i i ) . (14)
In practice, the parameter λ is first set to zero. Then, if the minimization algorithm [START_REF] Schultz | Model updating using constrained eigenstructure assignment[END_REF] does not converge, this value is increased until the algorithm converges. Let [S] with entries S i,j be the sensitivity matrix between the model parameters x j and the outputs z i . In the ideal case, after calculation of α 1 opt , . . . , α n opt , matrix [S] should be an identity matrix. If not, it means that the inputs are not identifiable separately and another function basis should be used. For instance, if a polynomial basis is chosen. Then the basis can be enriched with higher degree polynomials until [S] is an identity matrix. But it should be noted that increasing the size of the basis increases the size of matrices

[C], [C 1 ], . . . , [C n ].
Then if the size of the basis is too large, the calculation of α 1 , . . . , α n becomes computationally expensive.

Examples

Before introducing the model updating method in the next section, some examples are presented here in order to illustrate the construction of the transformed outputs. In the first example, x = (x 1 , x 2 , x 3 ) ∈ [1, 2] 3 and y(x) = (x 1 , x 2 + x 3 , 4 x 1 , x 3 , 3). The basis functions are such that f (y(x)) = y(x), i.e. z i (x) = α i T y(x) for i = 1, 2, 3. Therefore for each parameter x i , the transformed output will be a linear combination of the components of the output vector y(x). The previously described algorithm is run with λ = 5 × 10 -10 and 40, 000 Monte Carlo samples for the computation of the Sobol index matrices. The random variables X i have uniform distribution in [START_REF] Baruch | Optimisation procedure to correct stiffness and flexibility matrices using vibration data[END_REF][START_REF] Baruch | Modal data are insufficient for identification of both mass and stiffness matrices[END_REF] 3 . The initial values of α 1 , α 2 , α 3 are generated randomly. The following solutions are obtained: α 1 = (0.2, 0.0, 0.8, 0.0, 0.0), α 2 = (0.0, 3.4, 0.0, -3.4, 0.0), α 3 = (0.0, 0.0, 0.0, 3.4, 0.0), [START_REF] Sobol | Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates[END_REF] which correspond to the exact inverse reconstructions of x 1 , x 2 , x 3 respectively. For this case, the sensitivity matrix [S] is exactly identity.

In the second example, x = (x 1 , x 2 , x 3 ) ∈ [-0.4, 0.4] 3 and y(x) = (cos(x 1 + x 2 ), sin(x 2 ), x 1 x 2 , x 3 , 3). In a first case, the basis functions are, as the previous example, such that f (y(x)) = y(x), i.e.

As we can see in the first row of [S], the algorithm has failed to construct an output sensitive to x 1 only, which can be easily understood regarding the function y(x). For the second case, the function basis is enriched using using the squares of the model parameters in the combination, i.e. f (y(x)) = (y(x), y(x) 2 ). In this case, the new solutions are α 1 = (36.6, 0.0, 163.3, 0.0, 0.0, 68.1, 82.2, -5.8, 0.0, 0.0), α 2 = (0.0, 8.7, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0), α 3 = (0.0, 0.0, 0.0, 8.7, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0).

(

) 18
For this case, the matrix [S] is

[S] =   1.0 0.006 0.0 0.0 1.0 0.0 0.0 0.0 1.0   , (19)
therefore [S] is almost identity. The obtained result can be explained using quadratic approximations of cosine and sine functions which are present in the function y(x). For more complex functions more enrichments may be necessary.

Model updating

Once the vectors α 1 opt , . . . , α n opt are calculated, then the matrix [S] with entries S i,j can be estimated. Each entry S i,j represents how much the output z i (defined by vector α i ) is sensitive to the model parameter x j . Ideally, this matrix should be the identity matrix, meaning that for each model parameter x i , the function z i (constructed using the optimal value α i opt ) is exclusively sensitive to x i . In practice, this rarely happens (even if the basis is enriched) and, for each input x i , the discrepancy with respect to this ideal case can be quantified by the value Γ(α i opt ). If the value Γ(α i opt ) is sufficiently small then the model parameter x i can be identified independently of the other model parameters as the solution of the mean square problem,

x i opt = arg min x i ∈Cx i α i T f (y(x 1 , . . . , x i , . . . , x n )) -α i T f (y exp ) , (20)
in which α i T f (y exp ) is the transformed experimental output and where the parameters {x j , j = i} can have arbitrary values. We can then define a sequential algorithm for the identification of the model parameters: 1 -Rank the model parameters with respect to the values Γ(α i opt ). The indices of the ranked model parameters are denoted r 1 , . . . , r n . 2 -Keep only the n ǫ model parameters x r 1 , . . . , x rn ǫ for which Γ(α r i opt ) < ǫ, ǫ being a threshold to be defined.

3 -Initialize the model parameters vector x = x 0 . 4 -Identify x opt r 1 solving Eq.( 20) (using any univariate inverse problem solver) and replace x 0,r 1 by x opt r 1 . 5 -Identify iteratively x r 2 , . . . , x rn ǫ with the same method. Once all the model parameters x r 1 , . . . , x rn ǫ are identified, we can recalculate the vectors α r nǫ+1 , . . . , α rn of the remaining parameters (if any) and proceed with the same algorithm, the other model parameters being prescribed to their identified values. The algorithm has to be repeated until all the model parameters are identified or until all the Γ(α r i opt ) > ǫ which would mean that the remaining model parameters can't be identified separately.

The complete model updating method is described in Algorithm 1.

Algorithm 1 Model updating algorithm while n > 0 do Estimate [C], [C 1 ], . . . , [C n ]
Calculate α 1 opt , . . . , α n opt using Eq.( 14) Rank the model parameters such that Γ(α r1 opt ) < . . . < Γ(α rn opt ) Keep only the n ǫ model parameters for which Γ(α r i opt ) < ǫ Initialize the model parameters vector x opt for ℓ = 1, . . . , n ǫ do Identify x opt r ℓ solving Eq.( 20) (using any univariate inverse problem solver) Update the r th ℓ component of x opt end x ← (x nǫ+1 , . . . , x nǫ+1 ) n ← n -n ǫ end

Remarks i -Other measures of dependence

In the present paper, to construct outputs which depend only on one parameter, sensitivity indices have been used. Instead we could have used the mutual information (which is a measure of dependence between two random variables) between model parameters and the transformed output. Unfortunately this approach would require the estimation of probability distributions in high dimension for the calculation of the parameters α i opt and would not enable to separate α i as it was the case using sensitivity indices (see Eq.( 6)).

Another possibility would consist in weakening the independence using correlation between the system parameters and the output instead of the mutual information. This method is numerically very efficient but provides, as expected, good results only when the systems parameters and the outputs have approximatively Gaussian probability distributions (in which case correlation and dependence are equivalent).

ii -Computational cost

In the method presented here, the computational cost is mainly devoted to the estimation of matrices [C], [C 1 ], . . . , [C n ] using the Monte Carlo simulation method. This estimation is performed before the calculation of α 1 opt , . . . , α n opt . The other steps (calculation of α 1 opt , . . . , α n opt and identification of x opt ) are negligible.

iii -Modelling errors

If modelling errors are present due, for instance, to the simplicity of the model or measurement noise, then they can be taken into account using a probabilistic approach by introducing a random vector N at an advanced stage of the calculation. There are several ways to do this. For instance, an additive random vector N can be included directly in the output of the model:

ȳ(x) = y(x) + N . (21)
There are other ways to take into account uncertainties induced by modelling error. For instance, in structural dynamic, the nonparametric probabilistic approach [START_REF] Soize | Random matrix theory for modeling uncertainties in computational mechanics[END_REF] uses a random matrix approach to take into account these uncertainties. In this case, random vector N contains the entries of the random matrices. For all the cases, the matrix

[C j ] = [E X ∼j Cov X j (F | X ∼j ) ] is replaced by [ Cj ] = [E X ∼j Cov X j ,N (F | X ∼j ) ]. (22)
Then the rest of the method is kept unchanged. When introducing modelling errors, the variance of the output increases, making the search of a sensitivity matrix [S] equal to identity difficult. Therefore, the algorithm for calculating the vectors α i will automatically discard the components of the outputs that are highly impacted by modelling errors, making the updating method more robust to modelling errors.

Experimental application

In this section, the methodology presented in this paper is applied to the model updating of a 3-storey structure equipped with accelerators.

Test rig

The real structure is made up with 4 aluminium blocks connected each with the others with thin aluminium plates (see Fig. 1). The blocks have length 0.2 m, width 0.15 m and thickness 0.02 m. The thin plates have length 0.15 m, thickness 0.0012 m and width 0.015 m. The lowest block is clamped to the ground. The three other ones are equipped with accelerometers. The top block is excited with a hammer shock along x-direction. The measured Frequency Response Functions FRFs of the acceleration across x-direction are plotted in Fig. 2. As expected, we can see three peaks corresponding to the three global modes of the structure.

Model

The torsional and axial vibrations of the structure and the local vibration of the blocks and the plates are neglected. Therefore the structure is modelled as three masses connected by springs and dampers as shown in Fig. 3

ξ i = d i /(2 √ m i k i ), i = 1, 2, 3) are ξ 1 = ξ 2 = ξ 3 = 1.4 × 10 -3 .
The FRFs of the model (before updating) is compared to the experimental ones in Fig. 4. It can be seen in this figure large discrepancies in all the frequency band of analysis. To reduce this discrepancies, all the parameters will be updated in the next section.

Model updating

The model parameters that are used to update the model are: the three masses, the three stiffnesses and the three damping ratios. For each of these 9 parameters, a multiplicative updating coefficient is introduced. The prior output corresponds to the concatenation of the logarithm of the modulus of the FRFs shown in Fig. 4. A random noise is included in the model corresponding to independent and identically distributed centred Gaussian random variables with standard deviation 0.19. To construct the transformed outputs for the updating, we use λ = 1 × 10 -5 and 2, 000 Monte Carlo samples for the computation of the Sobol index matrices. The common logarithm (log 10 ) of the multiplicative updating parameters have uniform distributions in [-0.1, 0.1]. This parametrization enables the multiplicative parameters to vary between 0.1 and 10. The initial values of α i , i = 1, . . . , 9 are generated randomly. The so-obtained sensitivity matrix [S] is represented in Fig. 5. We can see in this figure that only the first six parameters (the three masses and the three stiffnesses) can be identified separately from the other parameters. We have tried to improved the sensitivity matrix by enriching the functional basis with quadratic and cubic functions. Nevertheless no noticeable improvement has been observed. Therefore the identification of the three other parameters (damping ratios) will require another iteration for which the results will be presented later. As described in Algorithm 1, the three masses and the three stiffnesses are then identified separately using the new outputs and a grid search for each parameter. The new values for these parameters are m 1 = 1.47 kg, m 2 = 1.52 kg, m 3 = 1.54 kg, k 1 = 1.41×10 4 N/m, k 2 = 1.30×10 4 N/m, k 3 = 1.41×10 4 N/m. The corresponding new FRFs are plotted in Fig. 6. Compared to the initial results plotted in Fig. 2, we can see that the updating of the masses and stiffnesses has drastically reduced the discrepancies between the experimental FRFs and the numerical one.

As explained before the damping ratios could not be identified in the previous iteration. Therefore, these parameters are calibrated in a new iteration by recalculating the vectors α i , with these parameters only. The value for the other parameters are fixed to the ones identified in the previous iteration. The sensitivity matrix [S] is represented in Fig. 7. The obtained matrix is not perfectly diagonal. But the error is sufficiently low to identify the damping ratios separately. The new values of the damping ratios are ξ 1 = 4.5 × 10 -3 , ξ 2 = 5.1×10 -4 and ξ 3 = 2.1×10 -3 . The corresponding new FRFs are plotted in Fig. 8. Compared to Fig. 6, we can see a slight improvement concerning the amplitude of the three resonance peaks.

Conclusions

In this paper, a new model updating method has been presented. This method enables the updating of each model parameter separately by constructing, for each parameter, a transform of the measured output which is sensitive to this parameter only. The sensitivity is measured by Sobol indices.

The separability of each parameter can be directly quantified using the obtained sensitivity matrix. Then each separable parameter can be updated using a univariate inverse problem solver. Then the procedure is restarted with the remaining non-separable parameters.

The method has been illustrated with a 3-storey structures for which experimental FRFs are measured. All the parameters of the model have been updated using two iterations of the algorithm. The response of the updated model is then very close to the experimental response. Such results would have been difficult to obtain using a manual updating.

In this work modelling errors are taken into account. Nevertheless other sources of uncertainties/variability may exist in the model. In this case a stochastic model should be constructed and should be updated using realizations of experimental outputs. The extension of the present work to the updating of such stochastic models will be addressed in future works.

  z i (x) = α i T y(x) for i = 1 . . . , 3. Again we use λ = 5 × 10 -10 and 40, 000 Monte Carlo samples for the computation of the Sobol index and the initial values of α 1 , α 2 , α 3 are generated randomly. The following solutions are obtained: α 1 = (33.0, 0.0, 84.2, 0.0, 0.0), α 2 = (0.0, 8.7, 0.0, 0.0, 0.0), α 3 = (0.0, 0.0, 0.0, 8.7, 0.0).

Figure 1 : 3 -

 13 Figure 1: 3-storey struture.

Figure 2 :

 2 Figure 2: Experimental FRF in acceleration for storey 1 (left) and storey 2 (right).

Figure 3 :

 3 Figure 3: Three masses model.

  . The masses of the plates are m 1 = m 2 = m 3 = 1.63 kg. The equivalent stiffnesses of the springs (representing the bending stiffness) are k 1 = k 2 = k 3 = 1.35 × 10 4 N/m. The damping ratio associated with the viscous dampers (

Figure 4 :

 4 Figure 4: Comparison between the experimental and numerical FRF in acceleration for storey 1 (left) and storey 2 (right).

Figure 5 :

 5 Figure 5: Sensitivity matrix [S].

Figure 6 :

 6 Figure 6: Comparison between the experimental and numerical FRF in acceleration for storey 1 (left) and storey 2 (right), after the first updating iteration (6 parameters only).

Figure 7 :

 7 Figure 7: Sensitivity matrix [S].

Figure 8 :

 8 Figure 8: Comparison between the experimental and numerical FRF in acceleration for storey 1 (left) and storey 2 (right), after the second updating iteration (3 damping ratios only).

Preprint submitted to Mechanical Systems and Signal ProcessingDecember 12, 2018


Hysteretic Bit/Rock Interaction Model to Analyze the Torsional Dynamics of a Drill String

Keywords: drill string nonlinear dynamics, bit/rock interaction model, torsional vibrations, stick-slip oscillations, stability map, hysteresis

The present paper proposes a novel hysteretic (non-reversible) bit/rock interaction model for the torsional dynamics of a drill string. Non-reversible means that the torque-on-bit depends not only on the bit speed, but also on the bit acceleration, producing a type of hysteretic cycle. The continuous drill string system is discretized by means of the finite element method and a reduced-order model is constructed using the normal modes of the associated conservative system. The parameters of the proposed hysteretic bit/rock interaction model is fitted with field data. The non-linear torsional vibration and the stability map of the drill string system are analyzed employing the proposed bit/rock interaction model and also a commonly used reversible model (without hysteresis). It turns out that the hysteretic model affects the stability region of the system.

Introduction

There are several papers available in the literature concerned with the drill string torsional dynamics and stick-slip oscillations [START_REF] Khulief | Vibration analysis of drillstrings with self-excited stick-slip oscillations[END_REF][START_REF] Navarro-López | Non-desired transitions and sliding-mode control of a multi-DOF mechanical system with stick-slip oscillations[END_REF][START_REF] Kreuzer | Controlling torsional vibrations of drill strings via decomposition of traveling waves[END_REF][START_REF] Ritto | Bayesian approach to identify the bit-rock interaction parameters of a drill-string dynamical model[END_REF][START_REF] Ritto | Validation of a drill string dynamical model and torsional stability[END_REF][START_REF] Mihajlovic | Friction-induced limit cycling in flexible rotor systems: An experimental drill-string set-up[END_REF][START_REF] Viguié | Using passive nonlinear targeted energy transfer to stabilize drill-string systems[END_REF]. A pure torsional model is sometimes enough to represent this kind of system. For instance, torsional models were applied successfully to represent test rigs that were constructed in [START_REF] Mihajlovic | Friction-induced limit cycling in flexible rotor systems: An experimental drill-string set-up[END_REF] to analyze the friction-induced limit cycling, and in [START_REF] Kreuzer | Controlling torsional vibrations of drill strings via decomposition of traveling waves[END_REF] to employ a control strategy.

There are many articles available in the literature concerned with friction laws, see for instance [START_REF] Wojewoda | Hysteretic effects of dry friction: modelling and experimental studies[END_REF][START_REF] Woodhouse | Are there reliable constitutive laws for dynamic friction?[END_REF]. The bit/rock interaction is very complex, but, in some situations, the relationship between the torque and the bit speed looks similar to a friction law. For instance, field measurements show (1) a torque weakening effect, and (2) hysteretic cycles [START_REF] Pavone | Application of high sampling rate downhole measurements for analysis and cure of stick-slip in drilling[END_REF][START_REF] Shen | The Origin and Mechanism of Severe Stick-Slip[END_REF]. Nevertheless, it should be emphasized that the bit/rock interaction is a much more complex mechanism.

Field data of a five kilometer drill string is analyzed in [START_REF] Ritto | Validation of a drill string dynamical model and torsional stability[END_REF], where again a pure torsional model presented satisfactory results reproducing field data, where torsional vibration was the dominant phenomenon observed. More generally, a coupled axial-lateral-torsional model should be applied [START_REF] Tucker | An integrated model for drill-string dynamics[END_REF][START_REF] Ritto | Bayesian approach to identify the bit-rock interaction parameters of a drill-string dynamical model[END_REF]. A full description model, including all dynamics, although possible, presents many difficulties due to lack of downhole data. During the drilling process there are many phenomena which are hard to measure, or simply not fully measured. Examples of these phenomenon include bit/rock interaction, fluid/rock interaction, proper well profile (inclination and azimuth), pipe/rock interaction, among others.

Figure 1 shows the field data [START_REF] Ritto | Validation of a drill string dynamical model and torsional stability[END_REF] that will support the proposed model presented in this paper. The downhole information used in this paper was acquired using a downhole mechanics measurement unit capable of providing both real-time measurement through mud telemetry and continuously recorded high-frequency data throughout the run. The sub, installed at the BHA above the bit, contains a suite of 19 sensors sampled at 10,000 Hz and downsampled and filtered prior to recording at 50 Hz. Following is a list of 50-Hz data recorded in this sub: triaxial accelerations; gyro rpm; magnetometer rpm; axial loading; torque; bending moment [START_REF] Shi | Integrated Efforts to Understand and Solve Challenges in 26-in Salt Drilling, Gulf of Mexico. SPE 180349-MS[END_REF].

The data show the dependence of the torque-on-bit with the bit speed; see Fig. 1. Closed to zero speed, the torque varies from about 6 kNm to 11 kNm, in a very steep straight line (close to 90 degrees). As the speed increases, the torque decreases. This type of torque-on-bit and bit speed 2 relation (weakening effect) has been observed experimentally before [START_REF] Pavone | Application of high sampling rate downhole measurements for analysis and cure of stick-slip in drilling[END_REF][START_REF] Shen | The Origin and Mechanism of Severe Stick-Slip[END_REF]. Some authors consider a non-linear function for the torque-on-bit versus bit speed to represent the bit rock/interaction model [START_REF] Khulief | Vibration analysis of drillstrings with self-excited stick-slip oscillations[END_REF][START_REF] Ritto | Bayesian approach to identify the bit-rock interaction parameters of a drill-string dynamical model[END_REF], while others apply a switching mechanism [START_REF] Leine | Stick-slip Whirl Interaction in Drillstring Dynamics[END_REF]. One can find coupled axial-torsional bit/rock interaction models, such as [START_REF] Richard | A simplified model to explore the root cause of stick-slip vibrations in drilling systems with drag bits[END_REF][START_REF] Tucker | Torsional vibration control and cosserat dynamics of a drill-rig assembly[END_REF].

Hysteretic cycles have been observed experimentally for the bit/rock interaction in [START_REF] Pavone | Application of high sampling rate downhole measurements for analysis and cure of stick-slip in drilling[END_REF][START_REF] Leine | Stick-slip Whirl Interaction in Drillstring Dynamics[END_REF]. Up to the authors knowledge the only hysteretic bit/rock interaction model found in the literature was proposed in [START_REF] Hong | Identification and control of stick?slip vibrations using Kalman estimator in oil-well drillstrings[END_REF]. The authors in [START_REF] Hong | Identification and control of stick?slip vibrations using Kalman estimator in oil-well drillstrings[END_REF] used the experimental results presented in [START_REF] Leine | Stick-slip Whirl Interaction in Drillstring Dynamics[END_REF], and applied their hysteretic model, which employs a switching mechanism, in the analysis of Proportional-Integral (PI) control strategy, aiming at mitigating stick-slip oscillations.

In the present paper a novel hysteretic (non-reversible) bit/rock interaction model is proposed based on the field data presented in [START_REF] Ritto | Validation of a drill string dynamical model and torsional stability[END_REF]. Nonreversible means that the torque-on-bit depends not only on the bit speed, but also on the bit acceleration, producing a type of hysteretic cycle. We call it hysteresis, even though when unloaded the torque goes back to zero. We depict the available field data, fit the model parameters with them, and argue that a hysteretic model would be appropriate for the bit/rock interaction. Then, the torsional vibration and stability map are analyzed.

This article is organized as follows. In Section 2 the torsional dynamical model is presented. The continuous system is discretized by means of the finite element method and a reduced-order model is constructed using the normal modes of the associated conservative system. The proposed bit/rock interaction model including hysteresis is also presented in Section 2. In Section 3 this interaction model is compared with field data. The numerical analysis are presented in Section 4 and, finally, the concluding remarks are made in Section 5.

Dynamic Model

Torsional model

The drill string system is basically composed by drill pipes (DP) and the bottomhole-assembly (BHA), as it is schematically represented in Fig. 2. DP are slender tubes that can reach kilometers, while BHA is composed by thicker tubes (drill collars) together with the measurement/formation logging equipment and a drill bit on the bottom, and its length usually reaches hundreds of meters. A vertical wellbore is considered, and only torsional vibrations are taken into account in the modelling. In this investigation, it is assumed that the reaction forces generated during the contact between the column and the wellbore and the column with the casing are small. It is also assumed that lateral and axial vibrations are small, compared to torsional dynamics. A constant speed Ω is imposed at the top, and a reaction torque appears due to the bit/rock interaction. In the field, a control strategy, which tries to maintain the speed constant, is applied to the top drive. A PID-control law is commonly used, but there are different torsional control strategies found in the literature [START_REF] Patil | A comparative review of modelling and controlling torsional vibrations and experimentation using laboratory setups[END_REF]. If the control is able to keep an almost constant speed at the top, our modeling is representative.

The angular rotation about the longitudinal axis, θ(x, t), where x is the longitudinal coordinate, is the solution of the following equation of motion:

ρI p (x) ∂ 2 θ(x, t) ∂t 2 -G ∂ ∂x I p (x) ∂θ(x, t) ∂x = T (x, t), (1)
where, I p is the cross sectional polar moment of inertia, ρ and G are the density and shear modulus of the material of the column, and T is the torque per unit length. The boundary conditions at the top (constant speed) are given by:

θ(0, t) = Ωt θ(0, t) = Ω , (2)
and the initial conditions are given by θ(x, 0) = 0, θ(x, 0) = Ω ,

Unlike [START_REF] Ritto | Non-linear dynamics of a drill-string with uncertain model of the bit-rock interaction[END_REF][START_REF] Ritto | Stochastic drill-string with uncertainty on the imposed speed and on the bit-rock parameters[END_REF], the present strategy solves the system considering its rotational displacements about a rotating frame. Let θ rel (x, t) be the relative angular rotation in the rotating frame associated with the imposed angle at the top. We then have

θ(x, t) = Ωt + θ rel (x, t). ( 4
)
The equation of motion is discretized by means of the finite element method, where linear shape functions are employed [START_REF] Inman | Engineering Vibration[END_REF]. Let u(t) be the vector of nodal coordinates related to the discretization of the of the relative angular rotation θ rel (x, t). Then the final discretized system, including a proportional damping model, is given by

[M]ü(t) + [D] u(t) + [K]u(t) = T( u(t), ü(t)), (5)
where [M] is the mass matrix, [D] is the damping matrix, [K] is the stiffness matrix, and T( u(t), ü(t)) is the torque vector.

The initial conditions are

u(0) = 0, u(0) = -Ω1. ( 6
)
where 1 is a vector with all entries equal to one. All the components of the generalized torque vector are zero except the one corresponding to the drill bit (last node, x = L). The non-linear torque applied to the bit is denoted by T bit ( θbit (t), θbit (t)) and will be described in the next section.

The normal modes of the conservative homogeneous system are used to construct a reduced-order model, considering the boundary conditions fixed at the top and free at the bottom. Few modes of the reduced-order basis give good approximations for the system analyzed, if the nonlinearity is weak. More modes are needed if the nonlinearity is severe. Other strategies, such as Karhunen-Loève Decomposition, or Proper Orthogonal Decomposition (POD) [START_REF] Kerschen | The method of proper orthogonal decomposition for dynamical characterization and order reduction of mechanical systems: An overview[END_REF][START_REF] Ritto | A new measure of efficiency for model reduction: Application to a vibroimpact system[END_REF], might be used to construct reduced-order basis for nonlinear systems.

The m first eigenvalues 0 < λ 1 ≤ λ 2 ≤ . . . ≤ λ m associated with elastic modes {ϕ 1 , ϕ 2 , . . . , ϕ m } are solutions of the generalized eigenvalue problem

[K] ϕ = λ[M] ϕ . (7)
The reduced-order model is obtained by projecting the full computational model on the subspace spanned by the m first elastic modes calculated using Eq. [START_REF] Mihajlovic | Friction-induced limit cycling in flexible rotor systems: An experimental drill-string set-up[END_REF]. Let [Φ] be n×m matrix whose columns are the m first elastic modes. We can then introduce the following approximation

u(t) = [Φ] q(t), (8)
in which q(t) is the vector of the m generalized coordinates which are solution of the reduced matrix equation

[ M ]q(t) + [ D] q(t) + [ K]q(t) = T( q(t), q(t)), (9)
with the initial conditions

q(0) = 0, q(0) = -Ω [ M ] -1 [Φ] T [M]1. (10)
In these equations,

[ M ] = [Φ] T [M] [Φ], [ D] = [Φ] T [D] [Φ] and [ K] = [Φ] T [K] [Φ]
are the m × m mass, damping and stiffness reduced-order matrices, and

T( q(t), q(t)) = [Φ] T T([Φ] q(t), [Φ]q(t))
is the reduced-order torque vector. The set of equations ( 8), ( 9) and ( 10) can be solved using commonly used integration schemes, such as the Euler scheme or the Runge-Kutta scheme, for instance. Damping ratios ζ i are selected directly to each i-th mode in the reduced damping matrix, which is diagonal and has entries 2ζ i ω ni , where i = 1, .., n, and ω ni is the i-th natural frequency of the system. In the simulations, ζ 1 = 0.095, ζ 2,3 = 0.030, and ζ i = 0.020 for i > 3.

Bit/Rock interaction model

The field data used in this study corresponds to the drilling operation of an ultradeepwater well. Due to a research collaboration between service company and operator, the BHA was equipped with a high-frequency measurement device. The sub is capable of measuring axial forces, torque, bending moments, angular speed, and tri-axial accelerations. The data set used in this study corresponds to the drilling of about 3 meters, where only the on-bottom data (actual drilling) was used.

The field data shown in Figure 1 have been smoothed (using a time-sliding window average) to remove the measurement noise and split into stick-slip cycles. Figure 3 shows six field data stick-slip cycles. As indicate the arrows of the first graphic of Fig. 3, when the bit accelerates, the path from above occurs, and when the bit speed decreases, the path from the bottom occurs. There is a fluctuation of the hysteretic cycles. Based on field experience, the possible cause of the variability of the cycles is due to the heterogeneity of the rock formation drilled.

These cycles indicate that (1) for each cycle, the torque value for positive acceleration and for negative acceleration are not the same, this is the hysteretic (non-reversible) phenomenon and (2) each observed cycle is different from one another. The present paper tackles the first point, i.e., it proposes a bit/rock interaction model including hysteresis and representing the mean bit/rock interaction. The construction of a stochastic bit/rock interaction model that takes into account the fluctuation of the hysteretic cycle is out of the scope of this paper and will be pursued in a future work. The mean cycle is obtained by applying a sliding window average, at each bit speed, for the upper phases (positive acceleration) and the lower phases (negative acceleration) separately. Figure 4 shows the field stick-slip cycles and their mean. The mean of all cycles (upper and lower) is in blue, and the mean cycle is shown in green. The upper green curve is the mean of the upper cycles and the lower green curve is the mean of the lower cycles. If we take the green line shown in Fig. 4, i.e., the field data mean hysteretic cycle, it is remarked that the thickness of this cycle decreases as the bit speed increases. At some point, the upper and lower curves collapse.

A commonly used non-linear reversible bit/rock interaction model is given by [START_REF] Khulief | Vibration analysis of drillstrings with self-excited stick-slip oscillations[END_REF][START_REF] Piovan | Continuous models for drill-strings of the oil industry: Analysis of approaches and discretization schemes[END_REF][START_REF] Ritto | Stochastic drill-string with uncertainty on the imposed speed and on the bit-rock parameters[END_REF]]

T bit ( θbit ) = µW r tanh(b 1 θbit ) + b 2 θbit 1 + b 3 θ2 bit , (11)
where T bit is the torque-on-bit in N.m, θbit is the bit speed in rad/s. The parameters b 1 , b 2 , and b 3 depend on the bit and rock properties. The weighton-bit is W, the friction coefficient is µ, and the bit radius is r. The weighton-bit, as expressed in the oil industry, is the amount of downward force exerted on the drill bit.

Figure 5 shows separately the two terms of the reversible model, Eq. [START_REF] Piovan | Continuous models for drill-strings of the oil industry: Analysis of approaches and discretization schemes[END_REF]. The first term (tanh, in blue) increases fast, and reaches a limit torque value as the speed increases. The second term (fraction term, in blue) is responsible for the peak in the bit/rock interaction model, and it dies out as the bit speed increases. The reversible model, Eq. ( 11) is the sum of these two terms (magenta curve).

This model considers a steep straight line, when the bit speed is close to zero. Usually this strategy produces good results. One can find in the literature other models that take into account explicitly stick and slip phases, such as in [START_REF] Leine | Stick-slip Whirl Interaction in Drillstring Dynamics[END_REF].

Once fitted, the bit/rock interaction model described by Eq. ( 11) yields an acceptable agreement with the average experimental plots (blue plots in Fig. 4). Nevertheless it suffers from two drawbacks: (1) it is not flexible enough to fit correctly the experiments for both low and high velocities; and (2) it cannot generate hysteresis effects in order to separate the forward and backward behaviours (green plots in Fig. 4).

In order to improve the model shown in Eq. ( 11) we suggest two modifications to address these two issues: (1) the power exponents of the bit speed in the fraction term can vary in order to give the model more flexibility; and (2) the fraction term is modulated by an acceleration (and speed)-dependent factor to obtain different amplitudes for the forward and the backward phases.

The new bit/rock interaction model, proposed in the present paper and which contains these improvements, is presented in Eq. [START_REF] Richard | A simplified model to explore the root cause of stick-slip vibrations in drilling systems with drag bits[END_REF].

T bit ( θbit , θbit ) = b 0 tanh(b 1 θbit ) + b 2 | θbit | b 4 sign( θbit ) 1 + b 3 | θbit | b 5 (1 + H( θbit , θbit )) , (12)
for θbit > 0. It is usually observed that the reaction torque is much lower when the bit speed is negative. This is because the bit/cutters are not symmetrical. Actually, we do not have field data to support any interaction model for θbit < 0, but, in the cases analyzed in the present paper, the values of the bit speed were always equal or greater than zero.

Note that if b 4 = 1, b 5 = 2 and H = 0 we retrieve the bit/rock interaction model found in Eq. [START_REF] Piovan | Continuous models for drill-strings of the oil industry: Analysis of approaches and discretization schemes[END_REF]. The parameters b 4 and b 5 , with 0 < b 4 < b 5 , depend on the bit and rock properties. The parameter b 0 [N.m] represents the weight-on-bit, the friction coefficient and the bit radius (for instance, this parameter might be written as b 0 = µW r).

The proposed model, Eq. ( 12), allows a better fit with the field data points. H is the hysteretic function: if it is equal to zero, the model is reversible, with no hysteretic cycles. The hysteretic model is only activated if H is different from zero, where H is presented in Eq. [START_REF] Ritto | Non-linear dynamics of a drill-string with uncertain model of the bit-rock interaction[END_REF].

H( θbit , θbit ) = β 1 tanh(β 2 θbit ) sign( θbit ) , (13)
which means that the variation of the bit speed (bit angular acceleration) should be taken into account, and the hysteretic cycle is limited (1 ± β 1 ).

Figure 6 shows the hysteretic term as a function of the bit acceleration. It increases fast, and reaches a limit value as the acceleration increases/decreases. In the present case β 1 = 14%, i.e., the hysteretic cycle is within plus or minus 14% of the reversible model. In addition, note that (1 -H) is multiplying only the second term of the bit/rock interaction model (the fraction term), such that as the bit speed increases the thickness of the hysteretic cycle decreases, as it is observed in most field data cycles shown in Fig. 3, and also in the field data mean cycle shown in Fig. 4. Now the model including hysteresis is considered, Eq.12 with H = 0. The fitted parameters of the hysteretic function are: β 1 = 14%, β 2 = 10.6, with appropriate units. Figure 9 shows the field data upper and lower mean cycles together with the fitted hysteretic bit/rock interaction model. The green lines are the same field data green lines as shown in Fig. 4. Again, there is a very good agreement between the proposed bit/rock interaction model and the field data. A larger hysteretic cycle is also considered, with β 1 = 56%. This parameter was chosen such that all field data cycles (black lines in in Fig. 4) fit in. Finally, Fig. 10 shows the three bit/rock interaction models that will be used for computations: (1) fitted reversible model (magenta line), ( 2) fitted model with hysteresis effects (continuous black line), and (3) model with large hysteresis effects (dashed black line).

The parameters used in the computational model are the ones that were fitted with the available field data:

b 0 = 3478, b 1 = 938, b 2 = 2.56, b 3 = 0.

Numerical results

In this section the results obtained for the three bit/rock interaction models are analyzed (see Fig. 10): ( 1) reversible (magenta), (2) hysteretic β 1 = 14% (black), and (3) hysteretic β 1 = 56% (black). The column is discretized in 100 finite elements. An adaptive Euler scheme is implemented to approximate the solution of the ODE. It is easy to implement and it is useful to compute the acceleration at each time step (which is necessary for the hysteretic H function). The time step used ∆t = 5 × 10 -4 s and each time step of the explicit integration scheme is computed as follows:

1. qi+1 = qi + ∆t [ M ] -1 ( T( qi , qi ) -[ D] qi -[ K]q i ) 2. q i+1 = q i + ∆t qi 3. qi+1 = ( qi+1 -qi )/∆t
The values of the parameters of the drill string system analyzed in the present paper are the following ones. Material properties ρ = 7800kg/m 3 , E = 220GPa, ν = 0.29, G = E/(2(1 + ν)), length of the drill pipe 4733.60 m, outer and inner diameter of the drill pipe 0.140 m and 0.119 m, length of the BHA 466.45 m, outer and inner diameter of the BHA 0.1607 m and 0.0727 m. The first five natural frequencies of the fixed-free torsional system, obtained using the stiffness and mass matrices of the finite element model, are [0.1316, 0.4199, 0.7388, 1.0719, 1.4113] Hz, and the damping rates related to the first five normal models are [0.095, 0.02, 0.02, 0.02, 0.02] (all the remaining damping rates are set to 0.02).

Figure 11 shows the bit speed time evolution when the nominal surface speed equals to 120 RPM and W equals to 245 kN. The torsional oscillations are larger when employing the hysteretic bit/rock interaction model. All of the responses show stick-slip oscillations, where the bit speed is zero and the system accumulates strain energy up to the point that the bit speed is released and reaches more than two times the nominal surface speed. The frequency of stick-slip is around 0.128 Hz for model 1 (reversible), 0.126 Hz for model 2 (small reversibility), and 0.112 Hz for model 3 (large hysteresis effects). Hence, the stick-slip frequency is a little lower than the first torsional natural frequency of the system, and it decreases as the hysteresis of the bit/rock interaction increases.

In order to construct a stability map, the nominal surface speed is analyzed from 50 to 160 RPM, and the weight-on-bit from 5 to 60 klf (22 to 267 kN). The bit/rock interaction model shown in Eq. 12 was calibrated considering 245 kN, and the model assumes that T bit is linear with respect to the weight-on-bit, Eqs. ( 12) and [START_REF] Ritto | Non-linear dynamics of a drill-string with uncertain model of the bit-rock interaction[END_REF]. Thus, a coefficient is used to multiply T bit such that different values of W can be simulated. For example, if the W is 200 kN, then T bit must be multiplied by 200/245 = 0.816. Let us define the stick-slip severity factor as SS = (ω max bit -ω min bit )/(2Ω). If there is no torsional oscillations, SS = 0. If there is stick-slip oscillations (ω min bit = 0) and the maximum bit speed (ω max bit ) is two times the nominal surface speed (Ω), then SS = 1. If SS is lower than 0.5 the system will be considered stable, otherwise it is considered unstable.

The stability map is constructed as follows. For each pair (W,Ω) the bit speed is computed and SS is recorded for the steady state response. Figure 12 shows the stability map when employing the bit/rock interaction models 1 (reversible) and 2 (small hysteresis effects). It can be noted that the maximum SS value is greater for the reversible bit/rock interaction model (model 2), but, at the same time, the stability region (dark blue) is a little bigger when employing this model. Now let us analyze Fig. 13, which shows the stability map when employing the bit/rock interaction model 3 (large hysteresis effects). This last chart presents a stability region (dark blue region) bigger than the other two charts.

This means, on one hand, that the hysteresis in the bit/rock interaction model favours the stability of the system. It seems that when the bit speed decreases, and T bit is in the lower curve of the hysteretic cycle (see Fig. 10), it allows the system to escape from instability. On the other hand, it can also be obverved that the instability region (orange-red region) in Fig. 13 presents values of SS higher than the charts in Fig. 12. Hence, at the same time that the hysteresis in the bit/rock interaction model favours the stability of the system for low values of weight-on-bit and high values of nominal surface speed, it favours the aggravation of the torsional oscillations of the system for high values of weight-on-bit and low values of nominal surface speeds. Since T bit passes through the upper curve of the hysteretic cycle (see Fig. 10) the amplitude of the stick-slip oscillations increases as well.

Concluding Remarks

The present paper proposes a bit/rock interaction model including a type of hysteresis (non-reversibility). The proposed model is based on field data and a hysteretic model is applied to represent the bit/rock interaction model of a drill string system.

Only torsional vibrations are analyzed for a torsion bar system discretized by means of the finite element method. A reduced-order model was constructed to speed up the computations. The torsional dynamics of the system is analyzed, and the stability map shows that the system including a bit/rock interaction model with hysteresis effects (1) favours the stability of the system for high nominal surface speeds and low weight-on-bit, and (2) induces higher stick-slip oscillations for low nominal surface speeds and high weight-on-bit. The reason for that has to do with the hysteretic cycle of the torque-on-bit versus the bit speed.

Since this paper presents the first results obtained employing this model, more investigations must be made with the proposed hysteretic bit/rock interaction model. In addition, the authors are now working on a stochastic model to take into account the fluctuations of the hysteretic cycles, which are present in the field data.

Figure 1 :

 1 Figure 1: Field data: (a) measured torque at the BHA, very close to the drill bit T exp bit , versus speed of the bit θbit , and (b) zoom image.

Figure 2 :

 2 Figure 2: General scheme of a drill string system.

Figure 3 :

 3 Figure 3: Downhole field data, torque [N.m] versus angular speed [RPM] at the bit: six stick-slip cycles available. The path from the top occurs when the bit accelerates, and the path from the bottom occurs when the bit speed decreases. The direction of the cycle is the same for all six cycles.

Figure 4 :

 4 Figure 4: Downhole field data: six stick-slip cycles and their mean. The mean of all cycles (upper and lower) is in blue, and the mean cycle is shown in green.

Figure 5 :

 5 Figure 5: Terms of the bit/rock interaction reversible model, Eq. (11): tanh term in green, fraction term in blue, and complete reversible model in magenta.

Figure 6 :

 6 Figure 6: (a) Hysteretic term H as a function of the bit acceleration and (b) zoom image.

Figure 7

 7 Figure7shows the hysteretic bit/rock interaction model for positive bit speeds. The arrows indicate the path of the cycle. The curve in blue happens when the bit acceleration is positive, and the curve in green happens when the bit accelaration is negative.

Figure 7 :

 7 Figure 7: Hysteretic bit/rock interaction model. Blue line when the bit acceleration is positive and green line when the bit accelaration is negative. The black straight vertical line indicates when the acceleration switches from positive to negative values.

3 .Figure 8 :

 38 Figure 8: Field data vs. bit/rock interaction reversible model, (Eq. 11): Fitted model in magenta and mean field data in blue.

Figure 9 :

 9 Figure 9: Field data vs. hysteretic bit/rock interaction model (Eq. 13): fitted model in black, and mean cycle of field data in green.

  38, b 4 = 0.78, b 5 = 1.1 and β 2 = 10.6. For model number one H = 0, for model number two β 1 = 14%, and for model number three β 1 = 56%.

Figure 10 :

 10 Figure 10: Bit/Rock interaction models (Eqs. 11 and 13): (1) reversible (magenta), (2) hysteretic β 1 = 14% (black), and (3) hysteretic β 1 = 56% (black).

Figure 11 :

 11 Figure 11: Bit speed for the three models (120 RPM, 245 kN): (1) reversible (magenta), (2) hysteretic β 1 = 14% (black), and (3) hysteretic β 1 = 56% (red).

Figure 12 :

 12 Figure 12: Stability map for models 1 (reversible) and 2 (hysteretic β 1 = 14%).

Figure 13 :

 13 Figure 13: Stability map for the model 3 (hysteretic β 1 = 56%).

Acknowledgments

The fourth author would like to acknowledgment the finantial support of the Brazilian agencies CNPq, CAPES, and FAPERJ. The authors would also like to acknowledge Schlumberger Oilfield Services for publishing this article.