Asymptotic Stabilization of the Oilwell Drillstring Torsional and Axial Vibrations

1. Introduction and Motivation

The problem of vibration quenching of the oilwell drillstring is extremely important from both technological and economic points of view. Oil extraction industry has at least 150-200 years hence it seems senseless to follow the entire literature on the subject. We shall restrict therefore to the last 2-3 decades, when the distributed parameter description was introduced (according to our knowledge). We met this description in [1], then the model occurred in [2], being afterwards adapted in the PhD thesis [3] and in the subsequent papers [4,5,6]. In the same line we can mention the surveys [7,8] and, with Particular reference, the monograph [9]. More recently we can mention the papers [10,11,12]

The aforementioned references deal with the distributed parameter drillstring dynamics and with the feedback control of vibration quenching. It is generally accepted that the (nonlinear) vibrations are induced by the interaction rock-drilling bit which is described by various nonlinear models [7,8,9]. On the other hand, there exist other models of vibration ignition, taking into account the dynamics of the driving motors while considering a lumped parameter dynamics for the drillstring. We can cite here the PhD thesis [13] and the paper [14].

While the idea of discovering hidden oscillations of [14] looks interesting and very useful, we shall focus here on the first class of the aforementioned problems i.e. vibration quenching for distributed parameter models. Our motivations consists in the necessity of completing the model with neglected terms allowing to obtain rigorous mathematical results for vibration quenching and in considering both torsional and axial vibrations. As mentioned in [10], the torsional vibrations are most important for the drillstring. While the model structure allows to consider the two types of vibrations separately, the mathematical aspects are slightly different.

Last but not least, the stability analysis, more precisely the asymptotic stability analysis, requires introduction of the neutral functional differential equations satisfying the well known assumption of J. K. Hale [15] about the asymptotic, even strong stability in the delays, of their difference operator. Throughout the paper it will appear how this assumption is fulfilled within the framework of the model deduction.

Therefore, what is left of the paper is organized as follows. Firstly, it is presented the dynamics of both torsional and axial vibrations as resulting from the application of the Hamilton variational principle [16,17,18], all basic notations being explained in the Appendix. Next, the steady states are defined and computed, the deviations from the state are defined. The systems in deviations are written and the energy-based Lyapunov functionals are associated to them. The properties of the models and of the derivatives of the Lyapunov functionals allow obtaining the uniform Lyapunov stability and global boundedness of the trajectories. These properties are obtained by synthesizing linear feedback controllers with dynamics, making use of the positiveness theory i.e. the Yakubovich-Kalman-Popov lemma [19].

In order to obtain asymptotic stability via the Barbashin-Krasovskii-LaSalle invariance principle, systems of coupled delay differential and difference equations (recognized to be of neutral type [15], p. 301) are associated to the initial models described by partial differential equations of hyperbolic type with derivative boundary conditions. The solutions of the two mathematical objects are shown to be in a one-to-one correspondence and, consequently, any mathematical result obtained for one of them is projected back on the other. Therefore asymptotic stability obtained for the systems of coupled delay differential and difference equations will give the same property for the initial models described by partial differential equations of hyperbolic type with derivative boundary conditions. The asymptotic stability results will be obtained first for each type of vibrations model separately. However the system of torsional vibrations introduces a perturbation in the system modeling the axial vibrations. Consequently the asymptotic stability is afterwards obtained for the coupled systems using a Lyapunov functional which is a linear combination of those used for each system describing the drillstring vibrations (torsional and axial). The paper ends with conclusions and enumeration of possible extensions of the obtained results.

2. The Model of the Coupled Torsional and Axial Vibrations of the Drillstring in $1D$ Distributed Parameters

The practical experience shows that the torsional vibrations are the dominant vibrations in drillstring dynamics. There exist however a number of notable references displaying both types of vibrations coupled at the level of the drilling bit, with the torsional vibrations acting on the axial ones but without the reverse [7,8,9]. Here we do not intend to compete with them, our aim being to deduce both dynamics from the Hamilton variational principle. For a general explanation of this principle, the reader is sent to such basic texts as [16,20,21,22]. In the following the Hamilton principle will be applied to mechanical systems subject to several external forces and torques. As it will appear in the sequel, the “list” of the external forces and torques is of utmost importance in order to obtain an adequate mathematical model.

2.1 Consider first the torsional - main - vibrations and construct the integral of the Hamilton principle. The rotational kinetic energy is given by

$${\displaystyle T_{\theta }\left(t\right)=\frac{1}{2}\left[J_m{\dot{\theta }}_m{\left(t\right)}^2+J_b{\theta }_t{(L,t)}^2+{\int }_0^L\rho \left(s\right)I_p\left(s\right){\theta }_t{(s,t)}^2\mathrm{d}s\right].}$$

(1)

The potential energy due to the elastic strain is given by

$${\displaystyle U_{\theta }\left(t\right)=\frac{1}{2}{\int }_0^{L}G\left(s\right)I_p\left(s\right){\nu }_s{(s,t)}^2=\frac{1}{2}{\int }_0^{L}G\left(s\right)I_p\left(s\right){\theta }_s{(s,t)}^2\mathrm{d}s.}$$

(2)

Obviously we considered the general case of a space inhomogeneous material i.e depending on $s\in [0,L]$ - the linear coordinate along the drillstring - see Appendix A with paper’s notations. The “list” of the generalized torques performing torsion and rotation, together with their expressions is also given in the Appendix.

We write down the work of the rotation torques

$$\begin{array}{c}{\displaystyle W_{m\theta }\left(t\right)=({\tau }_a\left(t\right)+{\tau }_0\left(t\right)+{\tau }_2\left(t\right)){\theta }_m\left(t\right)+({\tau }_1\left(t\right)+{\tau }_2\left(t\right))\theta (0,t)+}\hfill \\ {\displaystyle +({\tau }_0^{\prime }\left(t\right)+{\tau }_b\left(t\right))\theta (L,t)+{\int }_0^L{\tau }_d(s,t)\theta (s,t)\mathrm{d}s.}\hfill \end{array}$$

(3)

Consequently the Hamiltonian of the rotational dynamics is defined by

$${\displaystyle I_{\theta }(t_1,t_2)={\int }_{t_1}^{t_2}(T_{\theta }\left(t\right)-U_{\theta }\left(t\right)+W_{m\theta }\left(t\right))\mathrm{d}t,}$$

(4)

where $t_1$, $t_2$ are two arbitrary time moments.

2.2 Consider now the axial vibration dynamics. Define the kinetic energy

$${\displaystyle T_H\left(t\right)=\frac{1}{2}\left[m_0{\dot{z}}_H{\left(t\right)}^2+m_b{z}_t{(L,t)}^2+{\int }_0^L\rho \left(s\right)\Gamma \left(s\right)z_t{(s,t)}^2\mathrm{d}s\right]}$$

(5)

and the potential energy of the elastic strain by stretching and compression

$${\displaystyle U_H\left(t\right)=\frac{1}{2}{\int }_0^{L}E\left(s\right)\Gamma \left(s\right){\zeta }_s{(s,t)}^2\mathrm{d}s=\frac{1}{2}{\int }_0^{L}E\left(s\right)\Gamma \left(s\right)z_s{(s,t)}^2\mathrm{d}s.}$$

(6)

The list of the generalized forces performing penetration and deformation work on the drillstring (along the vertical) is given in the Appendix.

We now write down the work of the vertical forces as follows

$$\begin{array}{c}{\displaystyle W_H\left(t\right)=(f_a\left(t\right)+f_0\left(t\right)+f_2\left(t\right))z_H\left(t\right)+(f_1\left(t\right)+f_2\left(t\right))z(0,t)+}\hfill \\ {\displaystyle +(f_0^{\prime }\left(t\right)+f_b\left(t\right))z(L,t)+{\int }_0^L{f}_d(s,t)z(s,t)\mathrm{d}s}\hfill \end{array}$$

(7)

and the Hamiltonian of the axial dynamics

$${\displaystyle I_H(t_1,t_2)={\int }_{t_1}^{t_2}(T_H\left(t\right)-U_H\left(t\right)+W_H\left(t\right))\mathrm{d}t.}$$

(8)

The Hamiltonian of the overall drillstring vibration dynamics will then be

$${\displaystyle I(t_1,t_2)=I_{\theta }(t_1,t_2)+I_H(t_1,t_2).}$$

(9)

2.3 Consider now the Euler-Lagrange variations of the generalized coordinates as follows

$$\begin{array}{c}{\displaystyle \theta (s,t)=\overline{\theta }(s,t)+{\epsilon }_{\theta }\vartheta (s,t),{\theta }_m\left(t\right)={\overline{\theta }}_m\left(t\right)+{\epsilon }_{\theta }\vartheta \left(t\right)}\hfill \\ {\displaystyle z(s,t)=\overline{z}(s,t)+{\epsilon }_H\zeta (s,t),z_H\left(t\right)={\overline{z}}_H\left(t\right)+{\epsilon }_H{\zeta }_H\left(t\right),}\hfill \end{array}$$

(10)

where the bar variables account for a basic motion. Substituting (10) in the Hamiltonians the perturbed Hamiltonians are obtained

$${\displaystyle I^{\epsilon }(t_1,t_2)=I_{\theta }^{{\epsilon }_{\theta }}(t_1,t_2)+I_H^{{\epsilon }_H}(t_1,t_2)}$$

(11)

as follows

$$\begin{array}{c}{\displaystyle I_{\theta }^{{\epsilon }_{\theta }}(t_1,t_2)=\frac{1}{2}{\int }_{t_1}^{t_2}[J_m{({\dot{\overline{\theta }}}_m\left(t\right)+{\epsilon }_{\theta }{\dot{\vartheta }}_m\left(t\right))}^2+J_b{({\overline{\theta }}_t(L,t)+{\epsilon }_{\theta }{\vartheta }_t(L,t))}^2]\mathrm{d}t+}\hfill \\ {\displaystyle +\frac{1}{2}{\int }_{t_1}^{t_2}{\int }_0^L{I}_p\left(s\right)[\rho \left(s\right){({\overline{\theta }}_t(s,t)+{\epsilon }_{\theta }{\vartheta }_t(s,t))}^2-G\left(s\right){({\overline{\theta }}_s(s,t)+{\epsilon }_{\theta }{\vartheta }_s(s,t))}^2]\mathrm{d}s\mathrm{d}t+}\hfill \\ {\displaystyle +{\int }_{t_1}^{t_2}[({\tau }_a\left(t\right)+{\tau }_0\left(t\right)+{\tau }_2\left(t\right))({\overline{\theta }}_m\left(t\right)+{\epsilon }_{\theta }{\vartheta }_m\left(t\right))+({\tau }_1\left(t\right)+{\tau }_2\left(t\right))(\theta (0,t)+{\epsilon }_{\theta }\vartheta (0,t))]\mathrm{d}t+}\hfill \\ {\displaystyle +{\int }_{t_1}^{t_2}\left[({\tau }_0^{\prime }\left(t\right)+{\tau }_b\left(t\right))(\overline{\theta }(L,t)+{\epsilon }_{\theta }\vartheta (L,t))+{\int }_0^L{\tau }_d(s,t)(\overline{\theta }(s,t)+{\epsilon }_{\theta }\vartheta (s,t))\mathrm{d}s\right]\mathrm{d}t}\hfill \end{array}$$

(12)

and, analogously

$$\begin{array}{c}{\displaystyle I_H^{{\epsilon }_H}(t_1,t_2)=\frac{1}{2}{\int }_{t_1}^{t_2}[m_0{({\dot{\overline{z}}}_H\left(t\right)+{\epsilon }_H{\dot{\zeta }}_H\left(t\right))}^2+m_b{({\overline{z}}_t(L,t)+{\epsilon }_H{\zeta }_t(L,t))}^2]\mathrm{d}t+}\hfill \\ {\displaystyle +\frac{1}{2}{\int }_{t_1}^{t_2}{\int }_0^L\Gamma \left(s\right)[\rho \left(s\right){({\overline{z}}_t(s,t)+{\epsilon }_H{\zeta }_t(s,t))}^2-E\left(s\right){({\overline{z}}_s(s,t)+{\epsilon }_H{\zeta }_s(s,t))}^2]\mathrm{d}s\mathrm{d}t+}\hfill \\ {\displaystyle +{\int }_{t_1}^{t_2}[(f_a\left(t\right)+f_0\left(t\right)+f_2\left(t\right))({\overline{z}}_H\left(t\right)+{\epsilon }_H{\zeta }_H\left(t\right))+(f_1\left(t\right)+f_2\left(t\right))(\overline{z}(0,t)+{\epsilon }_H\zeta (0,t))]\mathrm{d}t+}\hfill \\ {\displaystyle +{\int }_{t_1}^{t_2}\left[(f_0^{\prime }\left(t\right)+f_b\left(t\right))(\overline{z}(L,t)+{\epsilon }_H\zeta (L,t))+{\int }_0^L{f}_d(s,t)(\overline{z}(s,t)+{\epsilon }_H\zeta (s,t))\mathrm{d}s\right]\mathrm{d}t.}\hfill \end{array}$$

(13)

Formulae (12) and (13) show that $I^{\epsilon }(t_1,t_2)$ is quadratic in both ${\epsilon }_{\theta }$ and ${\epsilon }_H$; however the two $"\epsilon "$-s never mix in the aforementioned expressions.

2.4 We shall recall now the Hamilton variational principle after [16], p.113:

The motion of the mechanical system occurs in such a way that the definite integral (11) becomes stationary for arbitrary possible variations of the configuration of the system, provided the initial and final configurations are given.

The stationarity conditions of the integral are given by

$$\begin{array}{c}{\displaystyle \frac{\partial }{\partial {\epsilon }_{\theta }}{I}^{\epsilon }(t_1,t_2){|}_{{\epsilon }_{\theta }=0}=\frac{\partial }{\partial {\epsilon }_{\theta }}{I}_{\theta }^{{\epsilon }_{\theta }}(t_1,t_2)=0}\hfill \\ {\displaystyle \frac{\partial }{\partial {\epsilon }_H}{I}^{\epsilon }(t_1,t_2){|}_{{\epsilon }_H=0}=\frac{\partial }{\partial {\epsilon }_H}{I}_H^{{\epsilon }_H}(t_1,t_2)=0.}\hfill \end{array}$$

(14)

It can be seen also from (12), (13) that the quadratic terms (in ${\epsilon }_{\theta }^2$ and ${\epsilon }_H^2$ are positive. Therefore the Hessian matrix of $I^{\epsilon }(t_1,t_2)$ is diagonal and positive definite. The stationary “point” defined by (14) is a minimum. Also “prescribed initial and final configurations” means that all variations are 0 for $t=t_1$, $t=t_2$:

$$\begin{array}{c}{\displaystyle {\vartheta }_m\left(t_1\right)=\vartheta (0,t_1)=\vartheta (L,t_1)=0;\vartheta (s,t_1)\equiv 0,\forall s\in [0,L]}\hfill \\ {\displaystyle {\vartheta }_m\left(t_2\right)=\vartheta (0,t_2)=\vartheta (L,t_2)=0;\vartheta (s,t_2)\equiv 0,\forall s\in [0,L]}\hfill \\ {\displaystyle {\zeta }_H\left(t_1\right)=\zeta (0,t_1)=\zeta (L,t_1)=0;\zeta (s,t_1)\equiv 0,\forall s\in [0,L]}\hfill \\ {\displaystyle {\zeta }_H\left(t_2\right)=\zeta (0,t_2)=\zeta (L,t_2)=0;\zeta (s,t_2)\equiv 0,\forall s\in [0,L].}\hfill \end{array}$$

(15)

The two integrals $I_{\theta }^{{\epsilon }_{\theta }}(t_1,t_2)$ and $I_H^{{\epsilon }_H}(t_1,t_2)$ can be discussed separately, as (14) shows. Consider first the integral $I_{\theta }^{{\epsilon }_{\theta }}(t_1,t_2)$ to find

$$\begin{array}{c}{\displaystyle \frac{\partial }{\partial {\epsilon }_{\theta }}{I}^{\epsilon }(t_1,t_2){|}_{{\epsilon }_{\theta }=0}={\int }_{t_1}^{t_2}[J_m{\dot{\overline{\theta }}}_m\left(t\right){\dot{\vartheta }}_m\left(t\right)+J_b{\overline{\theta }}_t(L,t){\vartheta }_t(L,t)]\mathrm{d}t+}\hfill \\ {\displaystyle +{\int }_{t_1}^{t_2}{\int }_0^L{I}_p\left(s\right)\left[\rho \left(s\right){\overline{\theta }}_t(s,t){\vartheta }_t(s,t)-G\left(s\right){\overline{\theta }}_s(s,t){\vartheta }_s(s,t)\right]\mathrm{d}s\mathrm{d}t+}\hfill \\ {\displaystyle +{\int }_{t_1}^{t_2}[({\tau }_a\left(t\right)+{\tau }_0\left(t\right)+{\tau }_2\left(t\right)){\vartheta }_m\left(t\right)+({\tau }_1\left(t\right)+{\tau }_2\left(t\right))\vartheta (0,t)]\mathrm{d}t+}\hfill \\ {\displaystyle +{\int }_{t_1}^{t_2}\left[({\tau }_0^{\prime }\left(t\right)+{\tau }_b\left(t\right))\vartheta (L,t)+{\int }_0^L{\tau }_d(s,t)\vartheta (s,t)\mathrm{d}s\right]\mathrm{d}t.}\hfill \end{array}$$

After the application of Fubini theorem and several integrations by parts, the following equations are obtained for any admissible motion

$$\begin{array}{c}{\displaystyle -J_m{\ddot{\theta }}_m\left(t\right)+{\tau }_a\left(t\right)+{\tau }_0\left(t\right)+{\tau }_2\left(t\right)=0}\hfill \\ {\displaystyle -J_b{\theta }_{tt}(L,t)+{\tau }_0^{\prime }\left(t\right)+{\tau }_b\left(t\right)-G\left(L\right)I_p\left(L\right){\theta }_s(L,t)=0}\hfill \\ {\displaystyle {\tau }_1\left(t\right)+{\tau }_2\left(t\right)+G\left(0\right)I_p\left(0\right){\theta }_s(0,t)=0}\hfill \\ {\displaystyle {\tau }_d(s,t)-\rho \left(s\right)I_p\left(s\right){\theta }_{tt}(s,t)+{\left(G\left(s\right)I_p\left(s\right){\theta }_s(s,t)\right)}_s=0.}\hfill \end{array}$$

(16)

Substituting the various expressions for the torques in (16), the following boundary value problem is obtained

$$\begin{array}{c}{\displaystyle \rho \left(s\right)I_p\left(s\right){\theta }_{tt}+c_{\theta }\left(s\right)I_p\left(s\right){\theta }_t-{\left(G\left(s\right)I_p\left(s\right){\theta }_s\right)}_s=0}\hfill \\ {\displaystyle J_m{\ddot{\theta }}_m+c_0{\dot{\theta }}_m+c_l{\theta }_t(0,t)-{\tau }_a\left(t\right)=0}\hfill \\ {\displaystyle c_l({\dot{\theta }}_m\left(t\right)-{\theta }_t(0,t))+G\left(0\right)I_p\left(0\right){\theta }_s(0,t)=0}\hfill \\ {\displaystyle J_b{\theta }_{tt}(L,t)+c_v{\theta }_t(L,t)+T_b\left({\theta }_t(L,t)\right)+G\left(L\right)I_p\left(L\right){\theta }_s(L,t)=0.}\hfill \end{array}$$

(17)

This is a boundary value problem for a hyperbolic second order partial differential equation. It is a non-standard boundary value problem since the boundary conditions contain derivatives and are coupled to certain ordinary differential equations. The partial differential equation represents the model of a vibrating string. This follows from the choice of the potential energy (2). Other equations, modeling Euler-Bernoulli, Kelvin-Voigt or Timoshenko beams can be obtained by choosing the corresponding forms of the potential energy [23,24,25].

Consider now the integral $I_H^{{\epsilon }_H}(t_1,t_2)$ whose structure is analogous to that of $I_{\theta }^{{\epsilon }_{\theta }}(t_1,t_2)$. Again, application of the Fubini theorem, integrations by parts and substitution of the various forces will give the following boundary value problem describing the dynamics of the axial vibrations

$$\begin{array}{c}{\displaystyle \rho \left(s\right)\Gamma \left(s\right)z_{tt}+c_H\left(s\right)\Gamma \left(s\right)z_t-{\left(E\left(s\right)\Gamma \left(s\right)z_s\right)}_s=0}\hfill \\ {\displaystyle m_0{\ddot{z}}_H+c_0^{\prime }{\dot{z}}_H+c_l^{\prime }{z}_t(0,t)-f_a\left(t\right)=0}\hfill \\ {\displaystyle c_l^{\prime }({\dot{z}}_H\left(t\right)-z_t(0,t))+E\left(0\right)\Gamma \left(0\right)z_s(0,t)=0}\hfill \\ {\displaystyle m_b{z}_{tt}(L,t)+c_v^{\prime }{z}_t(L,t)+E\left(L\right)\Gamma \left(L\right)z_s(L,t)-2{\left(\mu \kappa R_b\right)}^{-1}{T}_b\left({\theta }_t(L,t)\right)=0.}\hfill \end{array}$$

(18)

A simple examination of (17) and (18) shows that the equations of (17) are quasi-linear because of the term $T_b\left({\theta }_t(L,t)\right)$ which is nonlinear; they contain the input (control) torque signal ${\tau }_a\left(t\right)$ which serves to control the torsional vibrations by a feedback control device. The torsional vibrations dynamics is not affected by the axial dynamics and can be examined separately of it. Unlikely, examination of (18) shows that the axial vibrations dynamics system is linear, controlled by the input (control) force signal $f_a\left(t\right)$ but is also subject to an external perturbation signal $2{\left(\mu \kappa R_b\right)}^{-1}{T}_b\left({\theta }_t(L,t)\right)$ arising from the torsional vibrations dynamics. Consequently, the axial vibrations can also be analyzed separately but taking into account the effect of the perturbation.

2.5 To complete the discussion we shall give some details concerning the load torque at the drilling bit as they can be met in the basic references. One of the simplest form is given in [2]

$${\displaystyle T_b(z_t(L,t),{\theta }_t(L,t))=\left\{\begin{array}{c}{\displaystyle \zeta z_t(L,t){\mathrm{e}}^{-\alpha {\theta }_t(L,t)},{\theta }_t(L,t)>0}\hfill \\ {\displaystyle -\zeta z_t(L,t),{\theta }_t(L,t)<0}\hfill \end{array}\right.}$$

(19)

This expression introduces a coupling between torsional and axial dynamics; the author substitutes $z_t(L,t)$ by its steady state value, achieving the dynamics decoupling. In [7,8,9] the following expression is used, based on the models including Coulomb friction

$$\begin{array}{c}{\displaystyle T_b\left(z\right)=c_{b}z+W_{ob}{R}_b{\mu }_b\left(z\right)sgn\left(z\right);}\hfill \\ {\displaystyle {\mu }_b\left(z\right)={\mu }_{cb}+({\mu }_{sb}-{\mu }_{cb})exp\{-({\gamma }_b/{\nu }_f)\left|z\right|\}.}\hfill \end{array}$$

(20)

where $z={\theta }_t(L,t)$. This expression contains a finite jump at $z=0$

$${\displaystyle \underset{z\to 0_{+}}{lim}{T}_b\left(z\right)=W_{ob}{R}_b{\mu }_{sb};\underset{z\to 0_{-}}{lim}{T}_b\left(z\right)=-W_{ob}{R}_b{\mu }_{sb}.}$$

(21)

Little is known about initial boundary value problems for $1D$ hyperbolic partial differential equations with discontinuous right hand side and, in order to achieve a rigorous mathematical treatment, discontinouous expressions for $T_b(·)$ will be avoided. In those references and surveys previously cited [7,8,9] the discontinuous models of $T_b(·)$ - obtained from the use of various friction models (Coulomb, Karnopp) - were replaced by the following expression, considered to be a suitable approximation of them

$${\displaystyle T_b\left(z\right)=\frac{2pkz}{z^2+k^2}}$$

(22)

In fact, as its graphical representation shows - Figure 1a, the approximation is suitable around $z=0$ (the discontinuity point) mainly

Indeed, $T_b\left(0\right)=0$ and also

$${\displaystyle T_b^{\prime }\left(z\right)=\frac{2pk(k^2-z^2)}{{(z^2+k^2)}^2}.}$$

(23)

Smaller is k, closer to the discontinuity at $z=0$ is this approximation (p is a scale factor which adjusts the value of the extremum values). However, outside the interval $(-k,k)$, the slope of $T_b\left(z\right)$ changes sign (the graphic approaches slowly 0 asymptotically for $z\to \pm \infty$) and can induce instabilities. The underlying idea is that while $T_b\left(z\right)$ is subject to a sector (Lurie-type) condition

$${\displaystyle 0<\frac{T_b\left(z\right)}{z}<\frac{2p}{k},}$$

(24)

the condition of the global Lipschitz type that is

$${\displaystyle 0<\frac{T_b\left(z_1\right)-T_b\left(z_2\right)}{z_1-z_2}<L,}$$

(25)

which is necessary for the global asymptotic stability of an imposed steady state, is valid only on a bounded interval; consequently asymptotic stability holds only on bounded domains of the state space.

In [26,27] the following expression is used

$${\displaystyle T_b\left(z\right)=\frac{2T_{dyn}}{\pi }\left({\alpha }_{1}z{\mathrm{e}}^{-{\alpha }_2\left|z\right|}+arctan\left({\alpha }_{3}z\right)\right)}$$

(26)

with $z={\theta }_t(L,t)$. We give below some of its properties. Firstly we observe that $T_b\left(z\right)=-T_b(-z)$ hence $T_b\left(z\right)$ is an odd function: its analysis for $z>0$ is sufficient, also it suffices to consider $\Phi \left(z\right)=(\pi /\left(2T_{dyn}\right)T_b\left(z\right)$. Further

$${\displaystyle \Phi \left(0\right)=0,\underset{z\to \infty }{lim}\Phi \left(z\right)=\frac{\pi }{2}.}$$

Consider now the derivative ${\Phi }^{\prime }\left(z\right)$

$${\displaystyle {\Phi }^{\prime }\left(z\right)={\alpha }_1{\mathrm{e}}^{-{\alpha }_{2}z}-{\alpha }_1{\alpha }_{2}z{\mathrm{e}}^{-{\alpha }_{2}z}+\frac{{\alpha }_3}{1+{\left({\alpha }_{3}z\right)}^2}}$$

(27)

and ${\Phi }^{\prime }\left(0\right)={\alpha }_1+{\alpha }_3>0$ hence $\Phi \left(z\right)$ is increasing at least on a certain interval around $z=0$. The zeros of ${\Phi }^{\prime }\left(z\right)$ follow from the equation

$${\displaystyle {\alpha }_1[1+{\left({\alpha }_3\right)}^2](1-{\alpha }_{2}z)+{\alpha }_3{\mathrm{e}}^{{\alpha }_{2}z}=0}$$

(28)

and, taking into account the values of ${\alpha }_i$ given in [26,27], equation (28) has two positive roots. As shown there, $T_b\left(z\right)$ given by (26) looks as in Figure 1b, displaying a maximum and a minimum which are quite close, then increasing very slowly up to $\pi /2$. It follows that choosing e.g. $z_1$ such that $\left|\Phi \right(z_1\left)\right|<\pi /2$, the following will hold

$${\displaystyle {\displaystyle 0<\frac{\Phi \left(z_1\right)-\Phi \left(z_2\right)}{z_1-z_2}<{\alpha }_1+{\alpha }_3}}$$

(29)

(after some elementary straightforward manipulation). In fact (29) holds provided not both $z_1$, $z_2$ are located where the slope of $\Phi (·)$ is negative. In practice this condition is always fulfilled since the operating steady state point is defined by $\overline{\omega }$ - the steady state angular speed - within the interval $(-{\omega }_M,{\omega }_M)$ of the nonlinearity $T_b\left(\omega \right)$, ${\omega }_M$ corresponding to the maximal value of the nonlinear function $T_b\left(\omega \right)$. In the following we shall refer to (26) for $T_b\left(\omega \right)$.

3. Cyclic Variables. Steady States. Energy Identity.

3.1 A simple examination of (17) and (18) shows that ${\theta }_m\left(t\right)$, $\theta (s,t)$, $z_H\left(t\right)$, $z(s,t)$ enter in these equations only by their time derivatives i.e. they are cyclic variables; the equations give information about the derivatives only while the angular and linear positions are irrelevant (the assertion obviously follows from physical considerations). We can therefore introduce the new variables

$$\begin{array}{c}{\displaystyle {\omega }_m:={\dot{\theta }}_m,\omega (s,t):={\theta }_t(s,t),w(s,t):=G\left(s\right)I_p\left(s\right){\theta }_s(s,t)}\hfill \\ {\displaystyle v_H:={\dot{z}}_H,v(s,t):=z_t(s,t),w^H(s,t):=E\left(s\right)\Gamma \left(s\right)z_s(s,t)}\hfill \end{array}$$

(30)

to obtain the models without the cyclic variables

$$\begin{array}{c}{\displaystyle \rho \left(s\right)I_p\left(s\right){\omega }_t+c_{\theta }\left(s\right)I_p\left(s\right)\omega -w_s=0}\hfill \\ {\displaystyle w_t-G\left(s\right)I_p\left(s\right){\omega }_s=0}\hfill \\ {\displaystyle J_m{\dot{\omega }}_m+c_0{\omega }_m+c_l\omega (0,t)-{\tau }_a\left(t\right)=0}\hfill \\ {\displaystyle c_l({\omega }_m\left(t\right)-\omega (0,t))+w(0,t)=0}\hfill \\ {\displaystyle J_b{\omega }_t(L,t)+c_v\omega (L,t)+w(L,t)+T_b\left(\omega (L,t)\right)=0}\hfill \end{array}$$

(31)

and

$$\begin{array}{c}{\displaystyle \rho \left(s\right)\Gamma \left(s\right)v_t+c_H\left(s\right)\Gamma \left(s\right)v-w_s^H=0}\hfill \\ {\displaystyle w_t^H-E\left(s\right)\Gamma \left(s\right)v_s=0}\hfill \\ {\displaystyle m_0{\dot{v}}_H+c_0^{\prime }{v}_H+c_l^{\prime }v(0,t)-f_a\left(t\right)=0}\hfill \\ {\displaystyle c_l^{\prime }(v_H\left(t\right)-v(0,t))+w^H(0,t)=0}\hfill \\ {\displaystyle m_b{v}_t(L,t)+c_v^{\prime }v(L,t)+w^H(L,t)-\frac{2}{\mu \kappa R_b}{T}_b\left(\omega (L,t)\right)=0.}\hfill \end{array}$$

(32)

3.2 The equations of the steady states of (31) and (32) are obtained by letting the time derivatives to 0 thus obtaining

$$\begin{array}{c}{\displaystyle c_{\theta }\left(s\right)I_p\left(s\right)\overline{\omega }\left(s\right)-{\overline{w}}_s=0,G\left(s\right)I_p\left(s\right){\overline{\omega }}_s=0}\hfill \\ {\displaystyle c_0{\overline{\omega }}_m+c_l\overline{\omega }\left(0\right)-{\overline{\tau }}_a=0,c_l({\overline{\omega }}_m-\overline{\omega }\left(0\right))+\overline{w}\left(0\right)=0}\hfill \\ {\displaystyle c_v\overline{\omega }\left(L\right)+\overline{w}\left(L\right)+T_b\left(\overline{\omega }\left(L\right)\right)=0}\hfill \end{array}$$

(33)

and

$$\begin{array}{c}{\displaystyle c_H\left(s\right)\Gamma \left(s\right)\overline{v}\left(s\right)-{\overline{w}}_s^H\left(s\right)=0,E\left(s\right)\Gamma \left(s\right){\overline{v}}_s\left(s\right)=0}\hfill \\ {\displaystyle c_0^{\prime }{\overline{v}}_H+c_l^{\prime }\overline{v}\left(0\right)-{\overline{f}}_a=0,c^{\prime }l({\overline{v}}_H-\overline{v}\left(0\right))+{\overline{w}}^H\left(0\right)=0}\hfill \\ {\displaystyle c_v^{\prime }\overline{v}\left(L\right)+{\overline{w}}^H\left(L\right)=\frac{2}{\mu \kappa R_b}{T}_b\left(\overline{\omega }\left(L\right)\right).}\hfill \end{array}$$

(34)

These two systems are also rather similar. We shall start the analysis with system (33), obtaining firstly that $\overline{\omega }\left(s\right)\equiv \overline{\omega }=\mathrm{const}$. Next

$${\displaystyle \overline{w}\left(s\right)=\overline{w}\left(0\right)+\left({\int }_0^s{c}_{\theta }\left(\lambda \right)I_p\left(\lambda \right)\mathrm{d}\lambda \right)\overline{\omega }=c_l(\overline{\omega }-{\overline{\omega }}_m)+\left({\int }_0^s{c}_{\theta }\left(\lambda \right)I_p\left(\lambda \right)\mathrm{d}\lambda \right)\overline{\omega }.}$$

Therefore

$${\displaystyle \overline{w}\left(L\right)=c_l(\overline{\omega }-{\overline{\omega }}_m)+\left({\int }_0^L{c}_{\theta }\left(\lambda \right)I_p\left(\lambda \right)\mathrm{d}\lambda \right)\overline{\omega }=-c_v\overline{\omega }-T_b\left(\overline{\omega }\right)}$$

and $\overline{w}\left(s\right)$ is re-written as

$${\displaystyle \overline{w}\left(s\right)=-\left(c_v+{\int }_s^L{c}_{\theta }\left(\lambda \right)I_p\left(\lambda \right)\mathrm{d}\lambda \right)\overline{\omega }-T_b\left(\overline{\omega }\right).}$$

We have still to examine

$$\begin{array}{c}{\displaystyle c_0{\overline{\omega }}_m+c_l\overline{\omega }-{\overline{\tau }}_a=0}\hfill \\ {\displaystyle c_l(\overline{\omega }-{\overline{\omega }}_m)+c_v\overline{\omega }+\left({\int }_0^L{c}_{\theta }\left(\lambda \right)I_p\left(\lambda \right)\mathrm{d}\lambda \right)\overline{\omega }+T_b\left(\overline{\omega }\right)=0,}\hfill \end{array}$$

that is two nonlinear equations with three unknowns. Usually the angular velocity $\overline{\omega }$ of the drillstring is imposed from technological considerations. Then the system above becomes a linear one in the unknowns ${\overline{\omega }}_m$ and ${\overline{\tau }}_a$ thus allowing to obtain the steady state values for (31) from (33)

$$\begin{array}{c}{\displaystyle {\overline{\tau }}_a=\left(c_0+c_l+\frac{c_0{c}_v}{c_l}+\frac{c_0}{c_l}{\int }_0^L{c}_{\theta }\left(\lambda \right)I_p\left(\lambda \right)\mathrm{d}\lambda \right)\overline{\omega }+\frac{c_0}{c_l}{T}_b\left(\overline{\omega }\right)}\hfill \\ {\displaystyle {\overline{\omega }}_m=\left(1+\frac{c_v}{c_l}+\frac{1}{c_l}{\int }_0^L{c}_{\theta }\left(\lambda \right)I_p\left(\lambda \right)\mathrm{d}\lambda \right)\overline{\omega }+\frac{1}{c_l}{T}_b\left(\overline{\omega }\right)}\hfill \\ {\displaystyle \overline{w}\left(s\right)=-\left(c_v+{\int }_0^L{c}_{\theta }\left(\lambda \right)I_p\left(\lambda \right)\mathrm{d}\lambda \right)\overline{\omega }-T_b\left(\overline{\omega }\right),0\le s\le L.}\hfill \end{array}$$

(35)

We turn now to (34) where $T_b\left(\overline{\omega }\right)$ is a known steady state forcing term “arising” from “outside” (the system of torsional vibrations). We have $\overline{v}\left(s\right)\equiv \overline{v}=\mathrm{const}$. Here also, the advancement (penetrating) linear velocity $\overline{v}$ is imposed from technological considerations. Repeating the steps of the computation for the torsional vibrations steady state we obtain the steady state values for (32) from (34)

$$\begin{array}{c}{\displaystyle {\overline{f}}_a=\left(c_0^{\prime }+c_l^{\prime }+\frac{c_0^{\prime }{c}_v^{\prime }}{c_l^{\prime }}+\frac{c_0^{\prime }}{c_l^{\prime }}{\int }_0^L{c}_H\left(\lambda \right)\Gamma \left(\lambda \right)\mathrm{d}\lambda \right)\overline{v}-\frac{c_0^{\prime }}{c_l^{\prime }}\frac{2}{\mu \kappa R_b}{T}_b\left(\overline{\omega }\right)}\hfill \\ {\displaystyle {\overline{v}}_H=\left(1+\frac{c_v^{\prime }}{c_l^{\prime }}+\frac{1}{c_l^{\prime }}{\int }_0^L{c}_H\left(\lambda \right)\Gamma \left(\lambda \right)\mathrm{d}\lambda \right)\overline{v}-\frac{1}{c_l^{\prime }}\frac{2}{\mu \kappa R_b}{T}_b\left(\overline{\omega }\right)}\hfill \\ {\displaystyle {\overline{w}}^H\left(s\right)=-\left(c_v^{\prime }+{\int }_0^L{c}_H\left(\lambda \right)\Gamma \left(\lambda \right)\mathrm{d}\lambda \right)\overline{v}+\frac{2}{\mu \kappa R_b}{T}_b\left(\overline{\omega }\right),0\le s\le L}\hfill \end{array}$$

(36)

3.3 We shall consider now the energy identities for systems (31) and (32). We start with system (31). Multiply the first partial differential equation by $\omega (s.t)$, the second one by ${\left(G\left(s\right)I_p\left(s\right)\right)}^{-1}w(s,t)$, add the two resulting equalities and integrate the with respect to s from 0 to L to obtain the energy identity

$$\begin{array}{c}{\displaystyle \frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}{\int }_0^L\left[\rho \left(s\right)I_p\left(s\right)\omega {(s,t)}^2+\frac{1}{G\left(s\right)I_p\left(s\right)}w{(s,t)}^2\right]\mathrm{d}s+{\int }_0^L{c}_{\theta }\left(s\right)I_p\left(s\right)\omega {(s,t)}^2\mathrm{d}s-}\hfill \\ {\displaystyle -\omega (L,t)w(L,t)+\omega (0,t)w(0,t)\equiv 0.}\hfill \end{array}$$

(37)

In the same way we obtain the energy identity for (32)

$$\begin{array}{c}{\displaystyle \frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}{\int }_0^L\left[\rho \left(s\right)\Gamma \left(s\right)v{(s,t)}^2+\frac{1}{E\left(s\right)\Gamma \left(s\right)}{w}^H{(s,t)}^2\right]\mathrm{d}s+{\int }_0^L{c}_H\left(s\right)\Gamma \left(s\right)v{(s,t)}^2\mathrm{d}s-}\hfill \\ {\displaystyle -v(L,t)w^H(L,t)+v(0,t)w^H(0,t)\equiv 0.}\hfill \end{array}$$

(38)

The next step is to substitute in (37) $w(L,t)$ and $w(0,t)$ by taking into account the boundary conditions of (31) thus obtaining the identity

$$\begin{array}{c}{\displaystyle \frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}\left[J_m{\omega }_m{\left(t\right)}^2+J_b\omega {(L,t)}^2\right]+}\hfill \\ {\displaystyle +\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}{\int }_0^L\left[\rho \left(s\right)I_p\left(s\right)\omega {(s,t)}^2+\frac{1}{G\left(s\right)I_p\left(s\right)}w{(s,t)}^2\right]\mathrm{d}s+{\int }_0^L{c}_{\theta }\left(s\right)I_p\left(s\right)\omega {(s,t)}^2\mathrm{d}s+}\hfill \\ {\displaystyle +c_0{\omega }_m{\left(t\right)}^2+c_l\omega {(0,t)}^2+c_v\omega {(L,t)}^2+\omega (L,t)T_b\left(\omega (L,t)\right)-{\omega }_m\left(t\right){\tau }_a\left(t\right)\equiv 0.}\hfill \end{array}$$

(39)

The same procedure is used for (38) by taking into account the boundary conditions in (32)

$$\begin{array}{c}{\displaystyle \frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}\left[m_0{v}_H{\left(t\right)}^2+m_{b}v{(L,t)}^2\right]+}\hfill \\ {\displaystyle +\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}{\int }_0^L\left[\rho \left(s\right)\Gamma \left(s\right)v{(s,t)}^2+\frac{1}{E\left(s\right)\Gamma \left(s\right)}{w}^H{(s,t)}^2\right]\mathrm{d}s+{\int }_0^L{c}_H\left(s\right)\Gamma \left(s\right)v{(s,t)}^2\mathrm{d}s+}\hfill \\ {\displaystyle +c_0^{\prime }{v}_H{\left(t\right)}^2+c_l^{\prime }v{(0,t)}^2+c_v^{\prime }v{(L,t)}^2-v_H\left(t\right)f_a\left(t\right)-\frac{2}{\mu \kappa R_b}v(L,t)T_b\left(\omega (L,t)\right)\equiv 0.}\hfill \end{array}$$

(40)

The importance of the identities (39) and (40) will appear when applying the Lyapunov approach for controller synthesis and for proving asymptotic stability for the equilibrium describing the operating point of the drilling system.

4. The Systems in Deviations and Their Control. Stability

The problem which is tackled in this section is the design of the stabilizing feedback controller for the torsional and axial dynamics respectively as well as the achieved in this way stability properties.

4.1 Consider first the dynamics of the torsional vibrations (31) and its steady state (35). Introduce first the deviations from the steady state of the variables in (31)

$$\begin{array}{c}{\displaystyle {\xi }_m\left(t\right):={\omega }_m\left(t\right)-\overline{\omega },{\mu }_a\left(t\right):={\tau }_a\left(t\right)-\overline{\tau };}\hfill \\ {\displaystyle {\xi }_{\omega }(s,t):=\omega (s,t)-\overline{\omega },{\xi }_w(s,t):=w(s,t)-\overline{w}\left(s\right)}\hfill \end{array}$$

(41)

where $\overline{\omega }$ is given and ${\overline{\omega }}_m$, $\overline{w}\left(s\right)$ and $\overline{\tau }$ are obtained from (35). The new variables satisfy the system in deviations below

$$\begin{array}{c}{\displaystyle \rho \left(s\right)I_p\left(s\right){\left({\xi }_{\omega }\right)}_t+c_{\theta }\left(s\right)I_p\left(s\right){\xi }_{\omega }-{\left({\xi }_w\right)}_s=0}\hfill \\ {\displaystyle {\left({\xi }_w\right)}_t-G\left(s\right)I_p\left(s\right){\left({\xi }_{\omega }\right)}_s=0}\hfill \\ {\displaystyle J_m{\dot{\xi }}_m+c_0{\xi }_m+c_l{\xi }_{\omega }(0,t)-{\mu }_a\left(t\right)=0}\hfill \\ {\displaystyle c_l({\xi }_m\left(t\right)-{\xi }_{\omega }(0,t))+{\xi }_w(0,t)=0}\hfill \\ {\displaystyle J_b\frac{\mathrm{d}}{\mathrm{d}t}{\xi }_{\omega }(L,t)+c_v{\xi }_{\omega }(L,t)+{\xi }_w(L,t)+T_b({\xi }_{\omega }(L,t)+\overline{\omega })-T_b\left(\overline{\omega }\right)=0}\hfill \end{array}$$

(42)

Looking at (39) we can write down immediately the analogous identity associated to (42)

$$\begin{array}{c}{\displaystyle \frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}\left[J_m{\xi }_m{\left(t\right)}^2+J_b{\xi }_{\omega }{(L,t)}^2\right]+}\hfill \\ {\displaystyle +\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}{\int }_0^L\left[\rho \left(s\right)I_p\left(s\right){\xi }_{\omega }{(s,t)}^2+\frac{1}{G\left(s\right)I_p\left(s\right)}{\xi }_w{(s,t)}^2\right]\mathrm{d}s+}\hfill \\ {\displaystyle +c_0{\xi }_m{\left(t\right)}^2+c_l{\xi }_{\omega }{(0,t)}^2+c_v{\xi }_{\omega }{(L,t)}^2+{\xi }_{\omega }(L,t)[T_b({\xi }_{\omega }(L,t)+\overline{\omega })-T_b\left(\overline{\omega }\right)]+}\hfill \\ {\displaystyle +{\int }_0^L{c}_{\theta }\left(s\right)I_p\left(s\right){\xi }_{\omega }{(s,t)}^2\mathrm{d}s-{\xi }_m\left(t\right){\tau }_a\left(t\right)\equiv 0.}\hfill \end{array}$$

(43)

Identity (43) suggests a feedback controller with its output containing ${\xi }_m\left(t\right)$. On the other hand since the driving mechanism of the drillstring is located at $s=0$, the only possible input - technically speaking - is ${\xi }_m\left(t\right)$. We shall thus define the following linear controller

$${\displaystyle {\dot{x}}_c^{\omega }=A_c^{\omega }{x}_c^{\omega }+b_c^{\omega }{\xi }_m\left(t\right);{\nu }_c^{\omega }\left(t\right)={\left(f_c^{\omega }\right)}^T{x}_c^{\omega }+{\gamma }_c^{\omega }{\xi }_m\left(t\right)}$$

(44)

under the following basic assumptions on it: it has inherent stability i.e. $A_c^{\omega }$ is a Hurwitz matrix and the transfer function $H_{\omega }\left(\lambda \right)={\gamma }_c^{\omega }+{\left(f_c^{\omega }\right)}^T{(\lambda I-A_c^{\omega })}^{-1}{b}_c^{\omega }$ is an irreducible rational function - equivalently $(A_c^{\omega },b_c^{\omega })$ is a controllable pair and $({\left(f_c^{\omega }\right)}^T,A_c^{\omega })$ is an observable pair.

Taking ${\mu }_a\left(t\right)=-{\nu }_c^{\omega }\left(t\right)$ i.e. a negative feedback connection, the identity (43) becomes Looking at (39) we can write down immediately the analogous identity associated to (42)

$$\begin{array}{c}{\displaystyle \frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}\left[J_m{\xi }_m{\left(t\right)}^2+J_b{\xi }_{\omega }{(L,t)}^2\right]+}\hfill \\ {\displaystyle +\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}{\int }_0^L\left[\rho \left(s\right)I_p\left(s\right){\xi }_{\omega }{(s,t)}^2+\frac{1}{G\left(s\right)I_p\left(s\right)}{\xi }_w{(s,t)}^2\right]\mathrm{d}s+}\hfill \\ {\displaystyle +c_0{\xi }_m{\left(t\right)}^2+c_l{\xi }_{\omega }{(0,t)}^2+c_v{\xi }_{\omega }{(L,t)}^2+{\xi }_{\omega }(L,t)[T_b({\xi }_{\omega }(L,t)+\overline{\omega })-T_b\left(\overline{\omega }\right)]+}\hfill \\ {\displaystyle +{\int }_0^L{c}_{\theta }\left(s\right)I_p\left(s\right){\xi }_{\omega }{(s,t)}^2\mathrm{d}s+{\xi }_m\left(t\right)({\left(f_c^{\omega }\right)}^T{x}_c^{\omega }+{\gamma }_c^{\omega }{\xi }_m\left(t\right))\equiv 0.}\hfill \end{array}$$

(45)

We are now in position to define the following Lyapunov functional - written along the solutions of (42)-(44)

$$\begin{array}{c}{\displaystyle {\mathcal{V}}_{\omega }^{★}\left(t\right)={\left(x_c^{\omega }\left(t\right)\right)}^T{P}_{\omega }{x}_c^{\omega }\left(t\right)+\frac{1}{2}\left[J_m{\xi }_m{\left(t\right)}^2+J_b{\xi }_{\omega }{(L,t)}^2\right]+}\hfill \\ {\displaystyle +\frac{1}{2}{\int }_0^L\left[\rho \left(s\right)I_p\left(s\right){\xi }_{\omega }{(s,t)}^2+\frac{1}{G\left(s\right)I_p\left(s\right)}{\xi }_w{(s,t)}^2\right]\mathrm{d}s,}\hfill \end{array}$$

(46)

where $P_{\omega }$ is a symmetric matrix which will be determined in the following. Differentiating ${\mathcal{V}}_{\omega }^{★}\left(t\right)$ and taking into account the identity (45) we obtain

$$\begin{array}{c}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}{\mathcal{V}}_{\omega }^{★}\left(t\right)={\left(x_c^{\omega }\left(t\right)\right)}^T{P}_{\omega }(A_c^{\omega }{x}_c^{\omega }\left(t\right)+b_c^{\omega }{\xi }_m\left(t\right))+{(A_c^{\omega }{x}_c^{\omega }\left(t\right)+b_c^{\omega }{\xi }_m\left(t\right))}^T{P}_{\omega }{x}_c^{\omega }\left(t\right)-}\hfill \\ {\displaystyle -{\xi }_m\left(t\right)({\left(f_c^{\omega }\right)}^T{x}_c^{\omega }+{\gamma }_c^{\omega }{\xi }_m\left(t\right))-c_0{\xi }_m{\left(t\right)}^2-c_l{\xi }_{\omega }{(0,t)}^2-c_v{\xi }_{\omega }{(L,t)}^2-}\hfill \\ {\displaystyle -{\xi }_{\omega }(L,t)[T_b(\overline{\omega }+{\xi }_{\omega }(L,t))-T_b\left(\overline{\omega }\right)]-{\int }_0^L{c}_{\theta }\left(s\right)I_p\left(s\right){\xi }_{\omega }{(s,t)}^2\mathrm{d}s.}\hfill \end{array}$$

(47)

The right hand side of (47) contains obviously negative terms including the term in $T_b(·)$ - see Figure 1b. We have to focus on the first three terms in (47) - the ones in $x_c^{\omega }$ and ${\xi }_m$ containing also the symmetric matrix $P_{\omega }$ whose choice and properties are yet to be determined. We turn to Yakubovich-Kalman-Popov lemma [19,28] whose useful for our conditions version is given in Appendix 2. In our case $M=0$ and $\mathcal{G}(x,\mu )=-\mu ({\gamma }_c^{\omega }\mu +{\left(f_c^{\omega }\right)}^{T}x)$ and the lemma reads as follows

There exists a symmetric matrix $P_{\omega }=P_{\omega }^T$ and a vector $r_{\omega }$ of appropriate dimensions such that the identity

$$\begin{array}{c}{\displaystyle {\left(x_c^{\omega }\right)}^T(A_c^{\omega }{x}_c^{\omega }+b_c^{\omega }{\xi }_m)+{(A_c^{\omega }{x}_c^{\omega }+b_c^{\omega }{\xi }_m)}^T{P}_{\omega }{x}_c^{\omega }-{\xi }_m({\gamma }_c^{\omega }{\xi }_m+{\left(f_c^{\omega }\right)}^T{x}_c^{\omega })\equiv }\hfill \\ {\displaystyle \equiv -|\sqrt{{\gamma }_c^{\omega }}{\xi }_m+r_{\omega }^T{x}_c^{\omega }{|}^2}\hfill \end{array}$$

(48)

holds for any real vector $x_c^{\omega }$ and real scalar ${\xi }_m$ if and only if the following frequency domain inequality

$${\displaystyle \Re e{H}_{\omega }(ı\lambda )={\gamma }_c^{\omega }+\Re e{\left(f_c^{\omega }\right)}^T{(ı\lambda I-A_c^{\omega })}^{-1}{b}_c^{\omega }\ge 0}$$

(49)

holds for all frequencies $\lambda \in {\overline{\mathbb{R}}}_{+}$

Assuming that the controller is such that (49) holds, we take the symmetric matrix $P_{\omega }$ in (46) as the one prescribed by (48). Therefore the derivative (47) of the Lyapunov functional (47) is at least non-positive.

We next show that this will imply Lyapunov stability in the metrics induced by the Lyapunov functional itself of the zero equilibrium of the closed loop system in deviations, defined by (42),(44). Therefore we focus again on matrix $P_{\omega }$. From the Yakubovich-Kalman-Popov lemma it is known that (48) holds for any reals $x_c^{\omega }$ and ${\xi }_m$. Assume now that (48) holds, in particular, for some vector $x_c^{\omega }\left(t\right)$ and scalar ${\xi }_m\left(t\right)$ satisfying

$${\displaystyle {\dot{x}}_c^{\omega }\left(t\right)\equiv A_c^{\omega }{x}_c^{\omega }\left(t\right)+b_c^{\omega }{\xi }_m\left(t\right)}$$

(50)

even if they are not variables satisfying (42),(44). As a more particular case, assume ${\xi }_m\left(t\right)\equiv 0$ and integrate (48) between 0 and t to obtain

$${\displaystyle {\left(x_c^{\omega }\left(t\right)\right)}^T{P}_{\omega }{x}_c^{\omega }\left(t\right)-{\left(x_c^{\omega }\left(0\right)\right)}^T{P}_{\omega }{x}_c^{\omega }\left(0\right)=-{\left|r_{\omega }^T{x}_c^{\omega }\left(0\right)\right|}^2\le 0,}$$

hence

$${\displaystyle {\left(x_c^{\omega }\left(t\right)\right)}^T{P}_{\omega }{x}_c^{\omega }\left(t\right)\le {\left(x_c^{\omega }\left(0\right)\right)}^T{P}_{\omega }{x}_c^{\omega }\left(0\right).}$$

Since $A_c^{\omega }$ is a Hurwitz matrix, ${lim}_{t\to \infty }{\left(x_c^{\omega }\left(t\right)\right)}^T{P}_{\omega }{x}_c^{\omega }\left(t\right)=0$, ${\left(x_c^{\omega }\left(0\right)\right)}^T{P}_{\omega }{x}_c^{\omega }\left(0\right)\ge 0$, $P_{\omega }$ thus resulting nonnegative definite because $x_c^{\omega }\left(t\right)$ is arbitrary with its initial condition. From ${\mathcal{V}}^{★}\left(t\right)\le {\mathcal{V}}^{★}\left(0\right)$ and $P_{\omega }\ge 0$ we deduce that

$$\begin{array}{c}{\displaystyle \frac{1}{2}\left\{J_m{\xi }_m{\left(t\right)}^2+J_b{\xi }_{\omega }{(L,t)}^2+{\int }_0^L\left[\rho \left(s\right)I_p\left(s\right){\xi }_{\omega }{(s,t)}^2+\frac{1}{G\left(s\right)I_p\left(s\right)}{\xi }_w{(s,t)}^2\right]\mathrm{d}s\right\}\le }\hfill \\ {\displaystyle \le {\mathcal{V}}^{★}\left(0\right)={\left(x_c^{\omega }\left(0\right)\right)}^T{P}_{\omega }{x}_c^{\omega }\left(0\right){\left(x_c^{\omega }\left(0\right)\right)}^T{P}_{\omega }{x}_c^{\omega }\left(0\right)+\frac{1}{2}\left[J_m{\xi }_m{\left(0\right)}^2+J_b{\xi }_{\omega }{(L,0)}^2\right]+}\hfill \\ {\displaystyle +\frac{1}{2}{\int }_0^L\left[\rho \left(s\right)I_p\left(s\right){\xi }_{\omega }{(s,0)}^2+\frac{1}{G\left(s\right)I_p\left(s\right)}.{\xi }_w{(s,0)}^2\right]\mathrm{d}s.}\hfill \end{array}$$

(51)

On the other hand, let $\{x_c^{\omega }\left(t\right),{\xi }_m\left(t\right)\}$ be components of a solution of (42),(44). Since $A_c^{\omega }$ is a Hurwitz matrix, there exist ${\alpha }_0>0$, ${\gamma }_0>0$ such that ${\mathrm{e}}^{A_c^{\omega }t|}\le {\gamma }_0{\mathrm{e}}^{-{\alpha }_{0}t}$. Therefore a straightforward manipulation will give

$${\displaystyle |x_c^{\omega }\left(t\right)|\le {\gamma }_0\left\{|x_c^{\omega }\left(0\right)|+\frac{1}{{\alpha }_0}|b_c^{\omega }|\underset{t}{sup}\left|{\xi }_m\left(t\right)\right|\right\}\le {\gamma }_0\left|x_c^{\omega }\left(0\right)\right|+{\gamma }_1\sqrt{{\mathcal{V}}^{★}\left(0\right)}}$$

(52)

and (51),(52) give the Lyapunov stability of the zero equilibrium of (42)-(44), also global boundedness of the solutions. We have proven in fact

Theorem 1.

Consider the closed loop controlled system (42),(44) of the controlled torsional vibrations dynamics of the drillstring under the following assumptions: i) $A_c^{\omega }$ is a Hurwitz matrix; ii) the frequency domain inequality (49) holds; iii) the nonlinearity $T_b(·)$ defined by (29) satisfies

$${\displaystyle 0\le \frac{T_b\left(z_1\right)-T_b\left(z_2\right)}{z_1-z_2}\le \frac{2T_{dyn}}{\pi }({\alpha }_1+{\alpha }_2)}$$

(53)

Then the zero solution of (42),(44) is Lyapunov stable in the metrics induced by the Lyapunov functional (46), stability expressed by (51)-(52).

4.2 Consider now the dynamics of the axial vibrations. We shall proceed by analogy to the previous case of the torsional vibrations and pointing out the differences. Consider therefore the dynamics of the axial vibrations (32) and its steady state (36). Introduce the deviations from the steady state of the variables in (32)

$$\begin{array}{c}{\displaystyle {\xi }_H\left(t\right):=v_H\left(t\right)-{\overline{v}}_H,{\mu }_H\left(t\right):=f_a\left(t\right)-{\overline{f}}_a;}\hfill \\ {\displaystyle {\xi }_v(s,t):=v(s,t)-\overline{v},{\xi }_w^H(s,t):=w^H(s,t)-{\overline{w}}^H\left(s\right),}\hfill \end{array}$$

(54)

where $\overline{v}$ is given and ${\overline{v}}_H$, ${\overline{w}}^H\left(s\right)$ and ${\overline{f}}_a$ are obtained from (36). The new variables satisfy the system in deviations below

$$\begin{array}{c}{\displaystyle \rho \left(s\right)\Gamma \left(s\right){\left({\xi }_v\right)}_t+c_H\left(s\right)\Gamma \left(s\right){\xi }_v-{\left({\xi }_w^H\right)}_s=0}\hfill \\ {\displaystyle {\left({\xi }_w^H\right)}_t-E\left(s\right)\Gamma \left(s\right){\left({\xi }_v\right)}_t=0}\hfill \\ {\displaystyle m_0{\dot{\xi }}_H+c_0^{\prime }{\xi }_H+c_l^{\prime }{\xi }_v(0,t)-{\mu }_H\left(t\right)=0}\hfill \\ {\displaystyle c_l^{\prime }({\xi }_H\left(t\right)-{\xi }_v(0,t))+{\xi }_w^H(0,t)=0}\hfill \\ {\displaystyle m_b\frac{\mathrm{d}}{\mathrm{d}t}{\xi }_v(L,t)+c_v^{\prime }{\xi }_v(L,t)+{\xi }_w^H(L,t)-\frac{2}{\mu \kappa R_b}[T_b(\overline{\omega }+{\xi }_{\omega }(L,t))-T_b\left(\overline{\omega }\right)]=0.}\hfill \end{array}$$

(55)

Looking at (40) we can write the analogous identity associated to (55)

$$\begin{array}{c}{\displaystyle \frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}\left[m_0{\xi }_H{\left(t\right)}^2+m_b{\xi }_v{(L,t)}^2\right]+}\hfill \\ {\displaystyle +\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}{\int }_0^L\left[\rho \left(s\right)\Gamma \left(s\right){\xi }_v{(s,t)}^2+\frac{1}{E\left(s\right)\Gamma \left(s\right)}{\xi }_w^H{(s,t)}^2\right]\mathrm{d}s+}\hfill \\ {\displaystyle +c_0^{\prime }{\xi }_H{\left(t\right)}^2+c_l^{\prime }{\xi }_v{(0,t)}^2+c_v^{\prime }{\xi }_v{(L,t)}^2-\frac{2}{\mu \kappa R_b}{\xi }_v(L,t)[T_b({\xi }_{\omega }(L,t)+\overline{\omega })-T_b\left(\overline{\omega }\right)]+}\hfill \\ {\displaystyle +{\int }_0^L{c}_H\left(s\right)\Gamma \left(s\right){\xi }_v{(s,t)}^2\mathrm{d}s-{\xi }_H\left(t\right){\mu }_H\left(t\right)\equiv 0.}\hfill \end{array}$$

(56)

Identity (56) suggests a feedback controller with its output containing ${\xi }_H\left(t\right)$ - the deviation from the steady state. On the other hand, since the driving mechanism of the drillstring is located at $s=0$, the only possible input - technologically speaking - is also ${\xi }_H\left(t\right)$. By analogy to the case of the torsional vibrations, we shall define the following linear controller

$${\displaystyle {\dot{x}}_c^H=A_c^H{x}_c^H+b_c^H{\xi }_H\left(t\right);{\nu }_c^H\left(t\right)={\left(f_c^H\right)}^T{x}_c^H+{\gamma }_c^H{\xi }_H\left(t\right)}$$

(57)

under the same basic assumptions as in the previous case of the torsional vibrations: it has inherent stability i.e. $A_c^H$ is a Hurwitz matrix and the transfer function $H_H\left(\lambda \right)={\gamma }_c^H+{\left(f_c^H\right)}^T{(\lambda I-A_c^H)}^{-1}{b}_c^H$ is an irreducible rational function - equivalently $(A_c^H,b_c^H)$ is a controllable pair and $({\left(f_c^H\right)}^T,A_c^H)$ is an observable pair.

Now, a slight difference with respect to the previous case has to be taken into account: system (55) is subject to an external perturbation $2{\left(\mu \kappa R_b\right)}^{-1}[T_b(\overline{\omega }+{\xi }_{\omega }(L,t))-T_b\left(\overline{\omega }\right)]$ “arising” from the dynamics of the torsional vibrations and, as follows from Theorem 1, is bounded for $t\ge 0$. From now on we shall be interested in the inherent stability of (55) and will assume for a while that this external perturbation is identically zero.

We take ${\mu }_H\left(t\right)=-{\nu }_c^H\left(t\right)$ - a negative feedback connection and (56) becomes

$$\begin{array}{c}{\displaystyle \frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}\left[m_0{\xi }_H{\left(t\right)}^2+m_b{\xi }_v{(L,t)}^2\right]+}\hfill \\ {\displaystyle +\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}{\int }_0^L\left[\rho \left(s\right)\Gamma \left(s\right){\xi }_v{(s,t)}^2+\frac{1}{E\left(s\right)\Gamma \left(s\right)}{\xi }_w^H{(s,t)}^2\right]\mathrm{d}s+}\hfill \\ {\displaystyle +c_0^{\prime }{\xi }_H{\left(t\right)}^2+c_l^{\prime }{\xi }_v{(0,t)}^2+c_v^{\prime }{\xi }_v{(L,t)}^2+{\int }_0^L{c}_H\left(s\right)\Gamma \left(s\right){\xi }_v{(s,t)}^2\mathrm{d}s+}\hfill \\ {\displaystyle +{\xi }_H\left(t\right)({\left(f_c^H\right)}^T{x}_c^H\left(t\right)+{\gamma }_c^H{\xi }_H\left(t\right))\equiv 0}\hfill \end{array}$$

(58)

(the perturbation term is now identically zero). We can now define the Lyapunov functional - written along the solutions of (55),(57)

$$\begin{array}{c}{\displaystyle {\mathcal{V}}_H^{★}\left(t\right)={\left(x_c^H\left(t\right)\right)}^T{P}_H{x}_c^H\left(t\right)+\frac{1}{2}\left[m_0{\xi }_H{\left(t\right)}^2+m_b{\xi }_v{(L,t)}^2\right]+}\hfill \\ {\displaystyle +\frac{1}{2}{\int }_0^L\left[\rho \left(s\right)\Gamma \left(s\right){\xi }_v{(s,t)}^2+\frac{1}{E\left(s\right)\Gamma \left(s\right)}{\xi }_w^H{(s,t)}^2\right]\mathrm{d}s,}\hfill \end{array}$$

(59)

where $P_H$ is a symmetric matrix prescribed by the Yakubovich-Kalman-Popov lemma. Taking the derivative of ${\mathcal{V}}_H^{★}\left(t\right)$ combined with the identity (56) - without the perturbation generated term - it follows that

$$\begin{array}{c}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}{\mathcal{V}}_H^{★}\left(t\right)={\left(x_c^H\left(t\right)\right)}^T{P}_H(A_c^H{x}_c^H\left(t\right)+b_c^H{\xi }_m\left(t\right))+{(A_c^H{x}_c^H\left(t\right)+b_c^H{\xi }_H\left(t\right))}^T{P}_H{x}_c^H\left(t\right)-}\hfill \\ {\displaystyle -{\xi }_H\left(t\right)({\left(f_c^H\right)}^T{x}_c^H+{\gamma }_c^H{\xi }_H\left(t\right))-c_0^{\prime }{\xi }_H{\left(t\right)}^2-c_l^{\prime }{\xi }_v{(0,t)}^2-c_v^{\prime }{\xi }_v{(L,t)}^2-}\hfill \\ {\displaystyle -{\int }_0^L{c}_H\left(s\right)\Gamma \left(s\right){\xi }_v{(s,t)}^2\mathrm{d}s\equiv -{|\sqrt{{\gamma }_c^H}{\xi }_H\left(t\right)+r_H^T{c}_c^H\left(t\right)|}^2-}\hfill \\ {\displaystyle -c_0^{\prime }{\xi }_H{\left(t\right)}^2-c_l^{\prime }{\xi }_v{(0,t)}^2-c_v^{\prime }{\xi }_v{(L,t)}^2-{\int }_0^L{c}_H\left(s\right)\Gamma \left(s\right){\xi }_v{(s,t)}^2\mathrm{d}s\le 0.}\hfill \end{array}$$

(60)

The last identity follows from the application of Yakubovich-Kalman-Popov lemma under the positiveness frequency domain inequality

$${\displaystyle \Re e{H}_H(ı\mu )={\gamma }_c^H+\Re e{\left(f_c^H\right)}^T{(ı\mu I-A_c^H)}^{-1}{b}_c^H\ge 0}$$

(61)

for all frequencies $\mu \in {\mathbb{R}}_{+}$ - see (49) for comparison. FRom now on the development follows the line of the previous, torsional vibrations case. From the fact that $A_c^H$ is a Hurwitz matrix, we obtain that $P_H\ge 0$ and, as in the previous cae, it is obtained Lyapunov stability in the sense of the metrics induced by the Lyapunov functional itself. Without reproducing the formulae which are analogous to (51)-(52) we give below the result

Theorem 2.

Consider the closed loop control system (51),(57) of the controlled axial vibrations dynamics under the following assumptions: i) $A_c^H$ is a Hurwitz matrix; ii) the frequency domain inequality (61) holds; iii) the perturbation term $2{\left(\mu \kappa R_b\right)}^{-1}[T_b(\overline{\omega }+{\xi }_{\omega }(L,t))-T_b\left(\overline{\omega }\right)]$ is identically zero. Then the zero solution of (51),(57) is Lyapunov stable in the metrics induced by the Lyapunov functional (59).

5. Asymptotic Stability and Other Asymptotic Properties

The Lyapunov stability property has been obtained in the previous section by using the energy Lyapunov functional. As well known from the publications of the classics of the Lyapunov theory such as N. G. Četaev, the energy Lyapunov function(al) is a “weak” one which is easily constructed but with its derivative along the system’s solutions being only non-positive.(51),(57) - see e.g. [29]. This fact among others led to the Barbashin-Krasovskii-LaSalle invariance principle e.g. [30,31]. In the infinite dimensional case however, its proof and application requires pre-compactness of system’s bounded trajectories [32,33] - not easy to prove (notwithstanding, in the finite dimensional case the bounded trajectories are also compact).

For this reason we shall follow a pathe which goes back to Myshkis [34] and Cooke [35]. Our experience in this approach, summarized in [36] spreads over several decades. It consists in associating to the basic system a system of functional differential equations - in most cases of neutral case - for which the Barbashin-Krasovskii-LaSalle invariance principle holds. Due to the one-to-one correspondence between the solutions of the two aforementioned mathematical objects (the initial-boundary value problem with derivative boundary conditions and the associated system of functional differential equations with deviated argument) the asymptotic stability obtained for one of them (the system of functional differential equations) is automatically projected back on the other. In the following this approach will be applied to the dynamics of the drillstring vibrations.

5.1. Asymptotic stability of the equilibrium in the torsional vibrations case

We may call this aspect asymptotic vibration quenching in the torsional case. Consider the closed loop (with controller) system described by (42),(44) - with ${\mu }_a\left(t\right)=-{\nu }_c^{\omega }\left(t\right)$ - under the assumptions of zero distributed damping i.e.$c_{\theta }\left(s\right)\equiv 0$ and constant material parameters i.e. $I_p\left(s\right)\equiv I_p$, $\rho \left(s\right)\equiv \rho$, $G\left(s\right)\equiv G$. We shall thus have the following initial-boundary value problem with constant coefficients

$$\begin{array}{c}{\displaystyle \rho I_p{\left({\xi }_{\omega }\right)}_t-{\left({\xi }_w\right)}_s=0,{\left({\xi }_w\right)}_t-G{I}_p{\left({\xi }_{\omega }\right)}_s=0}\hfill \\ {\displaystyle J_m{\dot{\xi }}_m+(c_0+{\gamma }_c^{\omega }){\xi }_m+{\left(f_c^{\omega }\right)}^T{x}_c^{\omega }+c_l{\xi }_{\omega }(0,t)=0}\hfill \\ {\displaystyle {\dot{x}}_c^{\omega }=A_c^{\omega }{x}_c^{\omega }+b_c\omega {\xi }_m}\hfill \\ {\displaystyle c_l({\xi }_m\left(t\right)-{\xi }_{\omega }(0,t))+{\xi }_w(0,t)=0}\hfill \\ {\displaystyle J_b\frac{\mathrm{d}}{\mathrm{d}t}{\xi }_{\omega }(L,t)+c_v{\xi }_{\omega }(L,t)+{\xi }_w(L,t)+T_b(\overline{\omega }+{\xi }_{\omega }(L,t))-T_b\left(\overline{\omega }\right)=0.}\hfill \end{array}$$

(62)

We shall follow the procedure suggested in [35] and fully proven and described in [36]. Introduce firstly the Riemann invariants from

$$\begin{array}{c}{\displaystyle {\xi }_{\omega }(s,t)=\frac{1}{\sqrt{2}}(r^{+}(s,t)+r^{-}(s,t)/q)}\hfill \\ {\displaystyle {\xi }_w(s,t)=\frac{1}{\sqrt{2}}(-q{r}^{+}(s,t)+r^{-}(s,t)),q:=I_p\sqrt{\rho G}.}\hfill \end{array}$$

(63)

to associate to (62) the partial differential equations of the Riemann invariants

$${\displaystyle r_t^{+}+\lambda r_s^{+}=0,r_t^{-}-\lambda r_s^{-}=0;\lambda =\sqrt{G/\rho }.}$$

(64)

Substituting also (63) in the boundary conditions of (62), the transformed boundary conditions are associated to (64)

$$\begin{array}{c}{\displaystyle J_m{\dot{\xi }}_m+(c_0+{\gamma }_c^{\omega }){\xi }_m+{\left(f_c^{\omega }\right)}^T{x}_c^{\omega }+\frac{c_l}{q\sqrt{2}}(q{r}^{+}(0,t)+r^{-}(0,t))=0}\hfill \\ {\displaystyle {\dot{x}}_c^{\omega }=A_c^{\omega }{x}_c^{\omega }+b_c\omega {\xi }_m}\hfill \\ {\displaystyle q(q+c_l)r^{+}(0,t)-(q-c_l)r^{-}(0,t)-c_{l}q\sqrt{2}{\xi }_m\left(t\right)=0}\hfill \\ {\displaystyle J_b{\dot{\xi }}_b+c_v{\xi }_b+\frac{1}{\sqrt{2}}(r^{-}(L,t)-q{r}^{+}(L,t))+T_b(\overline{\omega }+{\xi }_b)-T_b\left(\overline{\omega }\right)=0}\hfill \\ {\displaystyle r^{-}(L,t)+q{r}^{+}(L,t)-q\sqrt{2}{\xi }_b\left(t\right)=0}\hfill \end{array}$$

(65)

(We denoted ${\xi }_{\omega }(L,t):={\xi }_b\left(t\right)$, the reason will appear in the sequel).

In the next step we shall consider the two families of the characteristic straight lines - the solutions of

$${\displaystyle \frac{\mathrm{d}t}{\mathrm{d}s}=\pm \frac{1}{\lambda }\Rightarrow t^{\pm }(\sigma ;s,t)=t\pm (\sigma -s)/\lambda }$$

(66)

and focusing on the two characteristics crossing some point $(s,t)\in [0,L]\times {\mathbb{R}}_{+}$. It is straightforward that the Riemann invariants are constant along the characteristics - $r^{+}$ along $t^{+}$ and $r^{-}$ along $t^{-}$. Therefore the following representation formulae are valid

$$\begin{array}{c}{\displaystyle r^{+}(s,t^{+}(s;s,t))\equiv r^{+}(s,t)=r^{+}(L,t^{+}(L;s,t))=r^{+}(L,t+(L-s)/\lambda )}\hfill \\ {\displaystyle r^{-}(s,t^{-}(s;s,t))\equiv r^{-}(s,t)=r^{-}(0,t^{-}(0;s,t))=r^{-}(0,t+s/\lambda ),}\hfill \end{array}$$

(67)

expressing the Riemann invariants as functions of their boundary values. Considering those characteristics which can be extended on the whole interval $[0,L]$, we obtain the following relations between the boundary values of the Riemann invariants

$${\displaystyle r^{+}(0,t)=r^{+}(L,t+L/\lambda ),r^{-}(L,t)=r^{-}(0,t+L/\lambda )}$$

(68)

and they emphasize why $T=L/\lambda$ is called the propagation time along the characteristics. Denoting

$${\displaystyle y^{+}\left(t\right):=r^{+}(L,t),y^{-}\left(t\right):=r^{-}(0,t);{\eta }^{\pm }\left(t\right):=y^{\pm }(t+T)}$$

(69)

and substituting ${\eta }^{\pm }\left(t\right)$ in (65) we obtain after some straightforward manipulation

$$\begin{array}{c}{\displaystyle J_m{\dot{\xi }}_m+\left(c_0+{\gamma }_c^{\omega }+\frac{c_l^2}{q+c_l}\right){\xi }_m+{\left(f_c^{\omega }\right)}^T{x}_c^{\omega }+\frac{c_l\sqrt{2}}{q+c_l}{\eta }^{-}(t-T)=0}\hfill \\ {\displaystyle {\dot{x}}_c^{\omega }=A_c^{\omega }{x}_c^{\omega }+b_c^{\omega }{\xi }_m}\hfill \\ {\displaystyle J_b{\dot{\xi }}_b+(q+c_v){\xi }_b-q\sqrt{2}{\eta }^{+}(t-T)+T_b(\overline{\omega }+{\xi }_b)-T_b\left(\overline{\omega }\right)=0}\hfill \\ {\displaystyle {\eta }^{+}\left(t\right)=\frac{c_l\sqrt{2}}{q+c_l}{\xi }_m\left(t\right)+\frac{q-c_l}{q(q+c_l)}{\eta }^{-}(t-T)}\hfill \\ {\displaystyle {\eta }^{-}\left(t\right)=-q{\eta }^{+}(t-T)+q\sqrt{2}{\xi }_b\left(t\right)}\hfill \end{array}$$

(70)

with $T:=L/\lambda$ the propagation time along the characteristics.

This system has the general form

$$\begin{array}{c}{\displaystyle \dot{x}\left(t\right)=Ax\left(t\right)+By(t-\tau )+f(x\left(t\right),y\left(t\right),y(t-\tau ))}\hfill \\ {\displaystyle y\left(t\right)=Cx\left(t\right)+Dy(t-\tau )+g\left(x\right(t),y(t),y(t-\tau \left)\right)}\hfill \end{array}$$

and was introduced in our early papers [37,38] and recognized as reducible to a system of neutral type in [15], p. 301. Besides the development of the second author of this paper - see [36] for a state of the art - this system has recently received a renewed approach e.g [39]. The solution of (70) can be constructed by steps on intervals $(kT,(k+1\left)T\right)$, $k\in \mathbb{Z}$ (i.e. forward and backward) provided initial conditions are given on $(-T,0)$. These initial conditions are obtained as follows: ${\xi }_m\left(0\right)$, $x_c^{\omega }\left(0\right)$, ${\xi }_b\left(0\right)$ “migrate” from (62) to (70) and ${\eta }^{\pm }\left(\theta \right)$ with $\theta \in [-T,0]$ follow from ${\xi }_{\omega }(s,0)$ and ${\xi }_w(s,0)$, $0\le s\le L$ via the converse of (63) - to obtain $r^{\pm }(s,0)$ - and from (67) by considering those characteristics which cannot be extended to the whole segment $[0,L]$ but cross it for some $s\in (0,L)$. In fact, from the constance of the Riemann invariants along the characteristics we shall have, based on (67), the following.From the first equality in (67), considering those characteristics crossing the abscissa $t=0$, it follows that the crossing takes place at $\widehat{\sigma }=s-\lambda t$, $0<s-\lambda t<L$. Therefore $r^{+}(L,t+(L-s)/\lambda )=r^{+}(s-\lambda t,0$ hence ${\eta }^{+}(t-s/\lambda )=r^{+}(s-\lambda t,0)$. In the same way, in the second equality of (67) we consider those characteristics crossing the abscissa $t=0$ at $\widehat{\sigma }=s+\lambda t$, $0<s+\lambda t<L$. Further we obtain ${\eta }^{-}(t+(s-L)/\lambda )=r^{-}(s+\lambda t,0)$. After a change of variables we obtain finally

$$\begin{array}{c}{\displaystyle {\eta }_0^{+}\left(\theta \right)=r^{+}(-\lambda \theta ,0)=\frac{1}{q\sqrt{2}}[q{\xi }_{\omega }(-\lambda \theta ,0)-{\xi }_w(-\lambda \theta ,0)]}\hfill \\ {\displaystyle {\eta }_0^{-}\left(\theta \right)=r^{-}(L+\lambda \theta ,0)=\frac{1}{\sqrt{2}}[q{\xi }_{\omega }(L+\lambda \theta ,0)+{\xi }_w(L+\lambda \theta ,0)]}\hfill \end{array}$$

(71)

It appears that starting from a classical solution of (62) with some initial conditions $\{{\xi }_m\left(0\right),x_c^{\omega }\left(0\right),{\xi }_b\left(0\right),{\xi }_{\omega }(s,0),{\xi }_w(s,0)\}$ we associated a solution of system (70) of functional differential equations with the initial conditions $\{{\xi }_m\left(0\right),x_c^{\omega }\left(0\right),{\xi }_b\left(0\right),{\eta }_0^{+}\left(\theta \right),{\eta }_0^{-}\left(\theta \right)\}$ with ${\eta }_0^{\pm }\left(\theta \right)$ obtained from (71). The converse is also true: considering a solution of (70) defined by some initial conditions $\{{\xi }_m\left(0\right),x_c^{\omega }\left(0\right),{\xi }_b\left(0\right),{\eta }_0^{\pm }\left(\theta \right)\}$ with ${\eta }_0^{\pm }\left(\theta \right)$ sufficiently smooth e.g. continuously differentiable, then $\{{\xi }_{\omega }(s,t),{\xi }_w(s,t),{\xi }_m\left(t\right),x_c^{\omega }\left(t\right),{\xi }_b\left(t\right)={\xi }_{\omega }(L,t)\}$ is a classical solution of (62) where $\{{\xi }_m\left(t\right),x_c^{\omega }\left(t\right),{\xi }_b\left(t\right)\}$ migrate from the solution of (70) and $\{{\xi }_{\omega }(s,t),{\xi }_w(s,t)\}$ are defined, as suggested by (67) and (69)

$${\displaystyle r^{+}(s,t)={\eta }^{+}(t-s/\lambda ),r^{-}(s,t)={\eta }^{-}(t+(s-L)/\lambda )}$$

(72)

and ${\xi }_{\omega }(s,t)$, ${\xi }_w(s,t)$ following from (63). The initial conditions of this solution follow by letting $t=0$ in the obtained expressions.

Summarizing we have obtained the following result

Theorem 3.

Consider a classical solution of (62) defined by some initial conditions $\{{\xi }_m\left(0\right),x_c^{\omega }\left(0\right),{\xi }_b\left(0\right),{\xi }_{\omega }(s,0),{\xi }_w(s,0)\}$, where we defined ${\xi }_b\left(t\right)={\xi }_{\omega }(L,t)$. Then $\{{\xi }_m\left(t\right),x_c^{\omega }\left(t\right),{\xi }_b\left(t\right),{\eta }^{\pm }\left(t\right)\}$ is a solution of system (70), where ${\eta }^{\pm }\left(t\right)$ are obtainedvia(63), (67), (69), is a solution of (70) with the initial conditions $\{\{{\xi }_m\left(0\right),x_c^{\omega }\left(0\right),{\xi }_b\left(0\right),{\eta }_0^{\pm }\left(\theta \right)\}$ where ${\eta }_0^{\pm }\left(\theta \right)$ are obtained from (71). Conversely, consider a solution of (70) defined by certain initial conditions $\{{\xi }_m\left(0\right),x_c^{\omega }\left(0\right),{\xi }_b\left(0\right),{\eta }_0^{\pm }(·)\}$ with ${\eta }_0^{\pm }(·)$ sufficiently smooth, then $\{{\xi }_{\omega }(s,t),{\xi }_w(s,t),{\xi }_m\left(t\right),x_c^{\omega }\left(t\right),{\xi }_b\left(t\right)={\xi }_{\omega }(L,t)\}$ is a classical solution of (62) with $\{{\xi }_{\omega }(s,t),{\xi }_w(s,t)\}$ defined from (72) and (63), the initial conditions being obtained by letting $t=0$ in the aforementioned solution of (62).

Now, Theorem 3 establishes a one-to-one correspondence between the solutions of two mathematical objects i.e. (62) and (70) and any result about one of these objects is projected back on the other.

In particular this will allow application of the Barbashin-Krasovskii-LaSalle invariance principle to system (70) to obtain asymptotic stability which will afterwards migrate to (62). To do this we shall associate to (70) a Lyapunov functional which is adapted from (46) - considered in the case of the constant parameters. In fact we have to express the integral in (46) using the representation formulae (72) and (63), recalling that $q=I_p\sqrt{\rho G}$ and $\lambda =\sqrt{G/\rho }$. The computations are quite straightforward and will give the following Lyapunov functional - written along the solutions of (70)

$$\begin{array}{c}{\displaystyle {\mathcal{V}}_{\omega }^{★}\left(t\right)=x_c^{\omega }{\left(t\right)}^T{P}_{\omega }{x}_c^{\omega }\left(t\right)+\frac{1}{2}\left[J_m{\xi }_m{\left(t\right)}^2+J_b{\xi }_b{\left(t\right)}^2\right]+}\hfill \\ {\displaystyle +\frac{1}{2}\left\{q{\int }_{-T}^0{\eta }^{+}{(t+\theta )}^2\mathrm{d}\theta +\frac{1}{q}{\int }_{-T}^0{\eta }^{-}{(t+\theta )}^2\mathrm{d}\theta \right\}.}\hfill \end{array}$$

(73)

Differentiating it along the solutions of (70) we obtain, based also on the re-writing of (47), relying on the representation formulae (72)

$$\begin{array}{c}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}{\mathcal{V}}_{\omega }^{★}\left(t\right)=-{|\sqrt{{\gamma }_c^{\omega }}{\xi }_m+r_{\omega }^T{x}_c^{\omega }|}^2-c_0{\xi }_m{\left(t\right)}^2-c_v{\xi }_b{\left(t\right)}^2-}\hfill \\ {\displaystyle -{\xi }_b\left(t\right)[T_b(\overline{\omega }+{\xi }_b\left(t\right))-T_b\left(\overline{\omega }\right)]-\frac{c_l}{2q^2}{(q{\eta }^{+}\left(t\right)+{\eta }^{-}(t-T))}^2\le 0.}\hfill \end{array}$$

(74)

This shows that the zero solution of (70) is Lyapunov stable in the sense of the metrics induced by the Lyapunov functional (73) itself. Based on Theorem 3 the Lyapunov stability in the sense of the Lyapunov functional (46) of system (42),(44) is also obtained. This last result however has been obtained directly in Section 4 in the more general case of the non-constant material parameters. But here we shall obtain also asymptotic stability by applying the Barbashin-Krasovskii-LaSalle invariance principle for neutral functional differential equations as it is stated in [15], p.293, Theorem 9.8.2.

We have to check firstly that the condition in Theorem 9.8.2 concerning the asymptotic stability of the difference operator is fulfilled. In the case of (70) this property is ensured by the location inside the unit disk of $\mathbb{C}$ of the eigenvalues of the matrix

$${\displaystyle D_{\omega }=\left(\begin{array}{ccc}0& & (q-c_l){\left(q(q+c_L)\right)}^{-1}\\ -q& & 0\end{array}\right)}$$

(75)

whose characteristic equation is

$${\displaystyle {\sigma }^2+\frac{q-c_l}{q+c_l}=0}$$

The property is true since $|(q-c_l){(q+c_l)}^{-1}|<1$. The Barbashin-Krasovskii-LaSalle invariance principle can thus be applied. In the set where $\mathrm{d}{\mathcal{V}}_{\omega }^{★}/\mathrm{d}t$ vanishes we have

$${\displaystyle r_{\omega }^T{x}_c^{\omega }\left(t\right)=0,{\xi }_m\left(t\right)=0,{\xi }_b\left(t\right)=0,q{\eta }^{+}\left(t\right)+{\eta }^{-}(t-T)=0}$$

(76)

It follows that in this set - the kernel of the derivative functional (74) - the dynamics of (70) is subject to the following system

$$\begin{array}{c}{\displaystyle {\dot{x}}_c^{\omega }=A_c^{\omega }{x}_c^{\omega },r_{\omega }^T{x}_c^{\omega }=0}\hfill \\ {\displaystyle J_m{\dot{\xi }}_m+{\left(f_c^{\omega }\right)}^T{x}_c^{\omega }+\frac{c_l\sqrt{2}}{q+c_l}{\eta }^{-}(t-T)=0}\hfill \\ {\displaystyle J_b{\dot{\xi }}_b-q\sqrt{2}{\eta }^{+}(t-T)=0}\hfill \\ {\displaystyle {\eta }^{+}\left(t\right)=\frac{q-c_l}{q(q+c_l)}{\eta }^{-}(t-T);{\eta }^{-}\left(t\right)=-q{\eta }^{+}(t-T).}\hfill \end{array}$$

(77)

From the first difference equation above and from (76) - the last equality - we obtain $2q{(q+c_l)}^{-1}{\eta }^{-}(t-T)=0$ hence ${\eta }^{-}\left(t\right)\equiv 0$ in any invariant set belonging to the kernel of the derivative. Therefore ${\eta }^{+}\left(t\right)\equiv 0$, ${\dot{\xi }}_m=0$, ${\dot{\xi }}_b=0$ hence ${\xi }_m\left(t\right)\equiv 0$, ${\dot{\xi }}_b\equiv 0$. Also since $x_c^{\omega }\left(t\right)\to 0$, we have ${\dot{x}}_c^{\omega }=0$ within the invariant set hence $x_c^{\omega }\left(t\right)\equiv 0$. Consequently the only invariant set included in the kernel of the derivative functional (74) is the equilibrium at 0. Application of the Barbashin- Krasovskii-LaSalle invariance principle will give, together with the Lyapunov stability already proven, the asymptotic stability of the zero equilibrium of the system of equations with deviated argument follows

$${\displaystyle \underset{t\to \infty }{lim}{x}_c^{\omega }\left(t\right)=0;\underset{t\to \infty }{lim}{\xi }_m\left(t\right)=\underset{t\to \infty }{lim}{\xi }_b\left(t\right)=\underset{t\to \infty }{lim}{\eta }^{\pm }\left(t\right)=0}$$

(78)

From the representation formulae (72) and from the definition of the Riemann invariants (63) it follows that

$${\displaystyle \underset{t\to \infty }{lim}{r}^{\pm }(s,t)=\underset{t\to \infty }{lim}{\xi }_{\omega }(s,t)=\underset{t\to \infty }{lim}{\xi }_w(s,t)=0,\forall s\in [0,L]\left(\mathrm{uniformly}\right)}$$

(79)

In fact we obtained asymptotic stability of the zero equilibrium of the system (62), the system of its Riemann invariants (65) and of the system of functional differential equations with deviated argument (70). The results of Section 4 and of SubSection 5.1 can be summarized as follows

Theorem 4.

Consider systems (62), (65), (70) under the following basic assumptions: i) all coefficients occurring in their equations are nonnegative; ii) the nonlinear function $T_b(·)$ is subject to a condition of (25) type; iii) the controller (44) has inherent stability i.e. matrix $A_c^{\omega }$ is a Hurwitz matrix and its transfer function $H_{\omega }\left(\lambda \right)$ is irreducible i.e. $(A_c^{\omega },b_c^{\omega })$ is controllable and $({\left(f_c^{\omega }\right)}^T,A_c^{\omega })$ is observable; iv) the transfer function $H_{\omega }\left(\lambda \right)$ satisfies the frequency domain inequality (49). Then the zero equilibrium of the aforementioned systems is globally asymptotically stable.

5.2. Asymptotic Properties of the Axial Vibrations Dynamics

5.2.1 Taking into account the similarities of the two kinds of dynamics, the exposition of the results in this subsection will be less extended than in the previous subsection. The assumptions on constant material parameters and zero distributed damping is valid here also. Consequently the basic closed loop system (i.e. with the synthesized controller) will result as follows

$$\begin{array}{c}{\displaystyle \rho \Gamma {\left({\xi }_v\right)}_t-{\left({\xi }_w^H\right)}_s=0,{\left({\xi }_w^H\right)}_t-E\Gamma {\left({\xi }_v\right)}_s=0}\hfill \\ {\displaystyle {\dot{x}}_c^H=A_c^H{x}_c^H+b_c^H{\xi }_H}\hfill \\ {\displaystyle m_0{\dot{\xi }}_H+(c_0^{\prime }+{\gamma }_c^H){\xi }_H+{\left(f_c^H\right)}^T{x}_c^H+c_l^{\prime }{\xi }_v(0,t)=0}\hfill \\ {\displaystyle c_l^{\prime }({\xi }_H\left(t\right)-{\xi }_v(0,t))+{\xi }_w^H(0,t)=0}\hfill \\ {\displaystyle m_b{\dot{\xi }}_b^H+c_v^{\prime }{\xi }_b^H+{\xi }_w^H(L,t)-2{\left(\mu \kappa R_b\right)}^{-1}[T_b(\overline{\omega }+{\xi }_{\omega }(L,t))-T_b\left(\overline{\omega }\right)]=0}\hfill \\ {\displaystyle {\xi }_v(L,t)={\xi }_b^H\left(t\right).}\hfill \end{array}$$

(80)

To this system we associate the Riemann invariants by

$$\begin{array}{c}{\displaystyle {\xi }_v(s,t)=\frac{1}{q_H\sqrt{2}}[q_H{r}_H^{+}(s,t)+r_H^{-}(s,t)]}\hfill \\ {\displaystyle {\xi }_w^H(s,t)=\frac{1}{\sqrt{2}}[-q_H{r}_H^{+}(s,t)+r_H^{-}(s,t)];q_H:=\Gamma \sqrt{\rho E}}\hfill \end{array}$$

(81)

and the system written in the Riemann invariants

$$\begin{array}{c}{\displaystyle {\left(r_H^{+}\right)}_t+{\lambda }_H{\left(r_h^{+}\right)}_s=0,{\left(r_H^{-}\right)}_t-\left({\lambda }_H\right){\left(r_H^{-}\right)}_s=0;{\lambda }_H=\sqrt{E/\rho }}\hfill \\ {\displaystyle {\dot{x}}_c^H=A_c^H{x}_c^H+b_c^H{\xi }_H}\hfill \\ {\displaystyle m_0{\dot{\xi }}_H+(c_0^{\prime }+{\gamma }_c^H){\xi }_H+{\left(f_c^H\right)}^T{x}_c^H+\frac{c_l^{\prime }}{q_H\sqrt{2}}[q_H{r}_H^{+}(0,t)+r_H^{-}(0,t)]=0}\hfill \\ {\displaystyle q_H(q_H+c_l^{\prime })r_H^{+}(0,t)-(q_H-c_l^{\prime })r_H^{-}(0,t)-c_l^{\prime }{q}_H\sqrt{2}{\xi }_H\left(t\right)=0}\hfill \\ {\displaystyle m_b{\dot{\xi }}_b^H+c_v^{\prime }{\xi }_b^H+\frac{1}{\sqrt{2}}[r_H^{-}(L,t)-q_H{r}_H^{+}(L,t)]-2{\left(\mu \kappa R_b\right)}^{-1}[T_b(\overline{\omega }+{\xi }_{\omega }(L,t))-T_b\left(\overline{\omega }\right)]=0}\hfill \\ {\displaystyle r_H^{-}(L,t)+q_H{r}_H^{+}(L,t)-q_H\sqrt{2}{\xi }_b^H\left(t\right)=0.}\hfill \end{array}$$

(82)

The two families of characteristic lines are solutions of the differential equations

$${\displaystyle \frac{\mathrm{d}t}{\mathrm{d}s}=\pm \frac{1}{{\lambda }_H}\Rightarrow t^{\pm }(\sigma ;s,t)=t\pm (\sigma -s)/{\lambda }_H.}$$

(83)

Focusing on the two characteristics crossing some point $(s,t)\in [0,L]\times {\mathbb{R}}_{+}$ and taking into account that the Riemann invariants are constant along the characteristics - $r^{+}$ along $t^{+}$ and $r^{-}$ along $t^{-}$ - the following representation formulae are obtained

$$\begin{array}{c}{\displaystyle r_H^{+}(s,t^{+}(s;s,t))=r_H^{+}(s,t)=r_H^{+}(L,t^{+}(L;s,t))=r_H^{+}(L,t+(L-s)/{\lambda }_H)}\hfill \\ {\displaystyle r_H^{-}(s,t^{-}(s;s,t))=r_H^{-}(s,t)=r_H^{-}(0,t^{-}(0;s,t))=r_H^{-}(0,t+s/{\lambda }_H).}\hfill \end{array}$$

(84)

giving the Riemann invariants as function of their boundary values. Considering those characteristics which can be extended on the whole interval $[0,L]$, we obtain the following relations between the boundary values of the Riemann invariants

$${\displaystyle r_H^{+}(0,t)=r_H^{+}(L,t+L/{\lambda }_H),r_H^{-}(L,t)=r_H^{-}(0,t+L/{\lambda }_H)}$$

(85)

and $T_H:=L/{\lambda }_H=L\sqrt{\rho /E}$ appears as the propagation time along the characteristics. Denoting

$${\displaystyle y_H^{+}\left(t\right):=r_H^{+}(L,t),y_H^{-}\left(t\right):=r_H^{-}(0,t);{\eta }_H^{\pm }\left(t\right):=y_H^{\pm }(t+T_H)}$$

(86)

and substituting ${\eta }_H^{\pm }\left(t\right)$ in (82), the following system of coupled delay differential and difference equations is obtained

$$\begin{array}{c}{\displaystyle {\dot{x}}_c^H=A_c^H{x}_c^H+b_c^H{\xi }_H}\hfill \\ {\displaystyle m_0{\dot{\xi }}_H+\left(c_0^{\prime }+{\gamma }_c^H+\frac{{\left(c_0^{\prime }\right)}^2}{q_H+c_l^{\prime }}\right){\xi }_H+{\left(f_c^H\right)}^T{x}_c+\frac{c_l^{\prime }\sqrt{2}}{q_H+c_l^{\prime }}{\eta }_H^{-}(t-T_H)=0}\hfill \\ {\displaystyle m_b{\dot{\xi }}_b^H+(q_H+c_v^{\prime }){\xi }_b^H-q_H\sqrt{2}{\eta }^{+}(t-T_H)-}\hfill \\ {\displaystyle -2{\left(\mu \kappa R_b\right)}^{-1}[T_b(\overline{\omega }+{\xi }_{\omega }(L,t))-T_b\left(\overline{\omega }\right)]=0}\hfill \\ {\displaystyle {\eta }_H^{+}\left(t\right)=\frac{c_l^{\prime }\sqrt{2}}{q_H+c_l^{\prime }}{\xi }_H\left(t\right)+\frac{q_H-c_l^{\prime }}{q_H(q_H+c_l^{\prime })}{\eta }_H^{-}(t-T_H)}\hfill \\ {\displaystyle {\eta }_H^{-}\left(t\right)=q_H\sqrt{2}{\xi }_H\left(t\right)-q_H{\eta }^{+}(t-T_H).}\hfill \end{array}$$

(87)

The solution of (87) can be constructed by steps based on the initial conditions $\{x_c^H\left(0\right),{\xi }_H\left(0\right),{\xi }_b^H\left(0\right);{\eta }_{H0}^{\pm }(·)\}$, where $\{x_c^H\left(0\right),{\xi }_H\left(0\right),{\xi }_b^H\left(0\right)\}$ migrate from the initial conditions of (80) and ${\eta }_{H0}^{\pm }\left(\theta \right)$, $-T_H\le \theta \le 0$, are constructed starting from $\{{\xi }_v(s,0),{\xi }_w^H(s,0)\}$ as follows

$$\begin{array}{c}{\displaystyle {\eta }_{H0}^{+}\left(\theta \right)=r_H^{+}(-{\lambda }_H\theta ,0)=\frac{1}{q_H\sqrt{2}}[q_H{\xi }_v(-{\lambda }_H\theta ,0)-{\xi }_w^H(-{\lambda }_H\theta ,0)]}\hfill \\ {\displaystyle {\eta }_{H0}^{-}\left(\theta \right)=r_H^{-}(L+{\lambda }_H\theta )=\frac{1}{\sqrt{2}}[q_H\sqrt{2}[q_H{\xi }_v(-{\lambda }_H\theta ,0)+{\xi }_w^H(-{\lambda }_H\theta ,0)],}\hfill \\ {\displaystyle -T_H\le \theta \le 0.}\hfill \end{array}$$

(88)

Conversely, starting from a sufficiently smooth solution of (87), defining

$${\displaystyle r_H^{+}(s,t)={\eta }_H^{+}(t-s/{\lambda }_H),r_H^{-}(s,t)={\eta }_H^{-}(t+(s-L)/{\lambda }_H)}$$

(89)

and ${\xi }_v(s,t)$, ${\xi }_w^H(s,t)$ from (81), then $\{x_c^H\left(t\right),{\xi }_H\left(t\right),{\xi }_b^H\left(t\right),{\xi }_v(s,t),{\xi }_w^H(s,t)\}$ is a classical solution of (80). The mathematical result is like Theorem 3 and reads

Theorem 5.

Consider a classical solution of (80) defined by some initial conditions ${\xi }_v(s,0)$, ${\xi }_w^H(s,0)$, ${\xi }_H\left(0\right)$, $x_c^H\left(0\right)$, ${\xi }_b^H\left(0\right)$, where we defined ${\xi }_b^H\left(t\right):={\xi }_w^H(L,t)$. Then $x_c^H\left(t\right)$, ${\xi }_H\left(t\right)$, ${\xi }_b^H\left(t\right)$, ${\eta }_H^{\pm }\left(t\right)$, where ${\eta }_H^{\pm }\left(t\right)$ are definedvia(81), (85), (86) is a solution of the system of functional differential equations (87) with the initial conditions ${\xi }_H\left(0\right)$, $x_c^H\left(0\right)$, ${\xi }_b^H\left(0\right)$, ${\eta }_{H0}^{\pm }\left(\theta \right)$, $-T_H\le \theta \le 0$, ${\eta }_{H0}^{\pm }\left(\theta \right)$ being obtained from (88). Conversely, consider a solution of (87) defined by certain initial conditions ${\xi }_H\left(0\right)$, $x_c^H\left(0\right)$, ${\xi }_b^H\left(0\right)$, ${\eta }_{H0}^{\pm }\left(\theta \right)$, $-T_H\le \theta \le 0$ with ${\eta }_{H0}^{\pm }(·)$ sufficiently smooth (e.g. continuously differentiable). Then ${\xi }_v(s,t)$, ${\xi }_w^H(s,t)$, ${\xi }_H\left(t\right)$, $x_c^H\left(t\right)$, ${\xi }_b^H\left(t\right)$ is a classical solution of (80) with ${\xi }_v(s,t)$, ${\xi }_w^H(s,t)$ defined from (89) and (81), the initial conditions resulting by letting $t=0$ in the aforementioned solution of (80). The forcing term $T_b(\overline{\omega }+{\xi }_{\omega }(L,t))-T_b\left(\overline{\omega }\right)$ is considered as given, being defined from the system describing the dynamics of the torsional vibrations.

As in the case of Theorem 3, Theorem 5 establishes a one-to-one correspondence between the solutions of two mathematical objects i.e. (80) and (87) and any result about one of these objects is projected back on the other. In particular this will allow application of the Barbashin- Krasovskii-LaSalle invariance principle to system (87) to obtain asymptotic stability which afterwards migrate to (80). Obviously the stability properties will be obtained for the autonomous systems i.e. by taking $T_b(\overline{\omega }+x{i}_{\omega }(L,t))-T_b\left(\overline{\omega }\right)\equiv 0$ in (80) and (87).

We shall thus associate to (87) a Lyapunov functional which is adapted from (59) in the case of constant parameters and by taking ${\xi }_v(L,t)={\xi }_b^H\left(t\right)$. Here also we have to express the integral in (59) using the representation formulae (81), (84), (86) and recalling that ${\lambda }_H=\sqrt{E/\rho }$, $q_H=\Gamma \sqrt{E\rho }$. The result will be

$$\begin{array}{c}{\displaystyle {\mathcal{V}}_H^{★}\left(t\right)=x_c^H{\left(t\right)}^T{P}_H{x}_c^H\left(t\right)+\frac{1}{2}\left[m_0{\xi }_H{\left(t\right)}^2+m_b{\xi }_b^H{\left(t\right)}^2\right]+}\hfill \\ {\displaystyle +\frac{1}{2}\left\{q_H{\int }_{-T_H}^0{\eta }_H^{+}{(t+\theta )}^2\mathrm{d}\theta +\frac{1}{q_H}{\int }_{-T_H}^0{\eta }^{-}{(t+\theta )}^2\mathrm{d}\theta \right\}.}\hfill \end{array}$$

(90)

Differentiating it along the solutions of (87) we obtain, based also on the re-writing of (60)

$$\begin{array}{c}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}{\mathcal{V}}_H^{★}\left(t\right)=-{|\sqrt{{\gamma }_c^H}{\xi }_H\left(t\right)+r_H^T{x}_c^H\left(t\right)|}^2-c_0^{\prime }{\xi }_H{\left(t\right)}^2-c_v^{\prime }{\xi }_b^H{\left(t\right)}^2-}\hfill \\ {\displaystyle -\frac{c_l^{\prime }}{2q_H^2}{[q_H{\eta }_H^{+}\left(t\right)+{\eta }_H^{-}(t-T_H)]}^2\le 0.}\hfill \end{array}$$

(91)

Inequality (91) shows that the zero solution of the autonomous (without forcing term) system (87) is Lyapunov stable in the sense of the metrics induced by the Lyapunov functional (90) itself. Based on Theorem 5, the Lyapunov stability in the metrics induced by the Lyapunov functional (59) for the system (80) is also obtained; this last result however was obtained directly in Section 4 in the most general case of the non-constant material parameters. However here we shall also obtain asymptotic stability by applying the Barbashin-Krasovskii-LaSalle invariance principle for neutral functional differential equations as it is stated in [15], p. 293, Theorem 9.8.2. here also, as in the case of the torsional vibrations, we have to check firstly that the condition in Theorem 9.8.2 concerning the asymptotic stability of the difference operator is fulfilled. In the case of system (87) this property is ensured by the location inside the unit disk of $\mathbb{C}$ of of the eigenvalues of the matrix

$${\displaystyle D_H=\left(\begin{array}{ccc}0& & (q_H-c_l^{\prime }){\left(q_H(q_H+c_L^{\prime })\right)}^{-1}\\ -q_H& & 0\end{array}\right)}$$

(92)

whose characteristic equation is

$${\displaystyle {\sigma }^2+\frac{q_H-c_l^{\prime }}{q_H+c_l^{\prime }}=0}$$

The property is true since $|(q_H-c_l^{\prime }){(q_H+c_l^{\prime })}^{-1}|<1$. The Barbashin-Krasovskii-LaSalle invariance principle can thus be applied. In the set where $\mathrm{d}{\mathcal{V}}_H^{★}/\mathrm{d}t$ vanishes we have

$${\displaystyle r_H^T{x}_c^H\left(t\right)=0,{\xi }_H\left(t\right)=0,{\xi }_b^H\left(t\right)=0,q_H{\eta }_H^{+}\left(t\right)+{\eta }_H^{-}(t-T_H)=0}$$

(93)

It follows that in this set - the kernel of the derivative functional (91) - the dynamics of (87) is subject to the following system

$$\begin{array}{c}{\displaystyle {\dot{x}}_c^H=A_c^H{x}_c^H,r_H^T{x}_c^H=0}\hfill \\ {\displaystyle m_0{\dot{\xi }}_H+{\left(f_c^H\right)}^T{x}_c^H+\frac{c_l^{\prime }\sqrt{2}}{q_H+c_l^{\prime }}{\eta }_H^{-}(t-T_H)=0}\hfill \\ {\displaystyle m_b{\dot{\xi }}_b^H-q_H\sqrt{2}{\eta }_H^{+}(t-T_H)=0}\hfill \\ {\displaystyle {\eta }_H^{+}\left(t\right)=\frac{q_H-c_l^{\prime }}{q_H(q_H+c_l^{\prime })}{\eta }_H^{-}(t-T_H);{\eta }_H^{-}\left(t\right)=-q_H{\eta }_H^{+}{(t-T)}_H.}\hfill \end{array}$$

(94)

From the first difference equation in (94) and from the last equality in (93) we obtain $2q_H{(q_H+c_l^{\prime })}^{-1}{\eta }_H^{-}(t-T_H)=0$ hence ${\eta }_H^{-}\left(t\right)\equiv 0$ in any invariant set included in the kernel of the derivative functional. Therefore ${\eta }_H^{+}\left(t\right)\equiv 0$, ${\dot{\xi }}_H\left(t\right)=0$, ${\dot{\xi }}_H^b\left(t\right)=0$ hence ${\xi }_H\left(t\right)\equiv 0$, ${\xi }_b^H\left(t\right)\equiv 0$. Also since $x_c^H\left(t\right)\to 0$, $A_c^H$ being a Hurwitz matrix, we have ${\dot{x}}_c^H\left(t\right)=0$, therefore $x_c^H\left(t\right)\equiv 0$.

Consequently the only invariant set included in the kernel of the derivative functional (91) is the equilibrium at 0. Application of the Barbashin-Krasovskii-LaSalle invariance principle will give, together with the Lyapunov stability already proven, asymptotic stability of the zero equilibrium of the autonomous system with deviated argument as follows

$${\displaystyle \underset{t\to \infty }{lim}{x}_c^H\left(t\right)=0;\underset{t\to \infty }{lim}{\xi }_H\left(t\right)=\underset{t\to \infty }{lim}{\xi }_b^H\left(t\right)=\underset{t\to \infty }{lim}{\eta }_H^{\pm }\left(t\right)=0}$$

(95)

From the representation formulae (89) and from the definition of the Riemann invariants (81) it follows that

$${\displaystyle \underset{t\to \infty }{lim}{r}_H^{\pm }(s,t)=\underset{t\to \infty }{lim}{\xi }_v(s,t)=\underset{t\to \infty }{lim}{\xi }_w^H(s,t)=0,\forall s\in [0,L]\left(\mathrm{uniformly}\right).}$$

(96)

In fact we obtained asymptotic stability of the zero equilibrium of the autonomous system (80), of the autonomous system (82) in the Riemann invariants and of the autonomous system (87) of functional differential equations with deviated argument. The results from Section 4 and of SubSection 5.2 can be summarized as follows

Theorem 6.

Consider the system (80) under the following basic assumptions: i) all the coefficients occurring in its equations are nonnegative; ii) the forcing exogeneous signal $T_b(\overline{\omega }+{\xi }_{\omega }(L,t))-T_b\left(\overline{\omega }\right)$ is identically zero; iii) the controller (57) has inherent stability i.e. matrix $A_c^H$ is a Hurwitz matrix, also its transfer function $H_H\left(\lambda \right)$ is irreducible i.e. $(A_c^H,b_c^H)$ is controllable and $({\left(f_c^H\right)}^T,A_c^H)$ is observable; iv) the transfer function $H_H\left(\lambda \right)$ satisfies the frequency domain inequality (61). Then the zero equilibrium of the aforementioned systems (80). (82) and (87) is globally asymptotically stable.

5.2.2 After the stability analysis of the autonomous system (80) and of the associated systems (82) and (87), it is still necessary to tackle the behavior under the perturbation $T_b(\overline{\omega }+{\xi }_{\omega }(L,t))-T_b\left(\overline{\omega }\right)$ which asymptotically approaches zero. The asymptotics of the perturbed systems has been studied for a long time. We can send the reader to such classical references as [30,31,40,41,42] dealing with stability with respect to persistent perturbations. Rather general and refined results on these topics can be found in the cycles of papers belonging to F. Brauer, A. Strauss and J. A. Yorke [43,44,45,46,47,48,49,50]. All aforementioned results deal with ordinary differential equations and strongly rely on the “good” properties of the associated Lyapunov functions.

Using these results here would probably require their extensions to functional differential equations or giving direct proofs based on the associated Lyapunov functionals. Fortunately, other approach, much simpler, is possible since the perturbing signal in (80) originates from system (62); these systems are connected in series - the signal $T_b(\overline{\omega }+{\xi }_{\omega }(L,t))-T_b\left(\overline{\omega }\right)$ is fed from (62) to (80) but no feedback exists between these systems. For these reasons we shall consider the composite system (62), (80) and the associated composite systems (65),(82) of the Riemann invariants and (70), (87) of the functional differential equations with deviated arguments. Taking into account the one-to-one correspondence between their solutions, we shall start the stability analysis by focusing on the last composite system of functional differential equations and associate the overall Lyapunov functional written along its solutions

$${\displaystyle {\mathcal{V}}^{★}\left(t\right)={\mathcal{V}}_{\omega }^{★}\left(t\right)+{\gamma }_0{\mathcal{V}}_H^{★}\left(t\right)}$$

(97)

where ${\mathcal{V}}_{\omega }^{★}\left(t\right)$ and ${\mathcal{V}}_H^{★}\left(t\right)$ are defined by (73) and (90) respectively; also $\gamma >0$. Taking the derivative of (98) along the solutions of (70) and (87) we shall obtain

$$\begin{array}{c}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}{\mathcal{V}}^{★}\left(t\right)=-{|\sqrt{{\gamma }_c^{\omega }}{\xi }_m\left(t\right)+r_{\omega }^T{x}_c^{\omega }\left(t\right)|}^2-c_0{\xi }_m{\left(t\right)}^2-c_v{\xi }_b{\left(t\right)}^2-}\hfill \\ {\displaystyle -{\xi }_b\left(t\right)[T_b(\overline{\omega }+{\xi }_b\left(t\right))-T_b\left(\overline{\omega }\right)]-\frac{c_l}{2q^2}{(q{\eta }^{+}\left(t\right)+{\eta }^{-}(t-T))}^2-}\hfill \\ {\displaystyle -{\gamma }_0\left\{|\sqrt{{\gamma }_c^H}{\xi }_H\left(t\right)+r_H^T{x}_c^H{\left(t\right)|}^2+c_0^{\prime }{\xi }_H{\left(t\right)}^2+c_v^{\prime }{\xi }_b^H{\left(t\right)}^2+\right.}\hfill \\ {\displaystyle +\left.\frac{c_l^{\prime }}{2q_H^2}{[q_H{\eta }_H^{+}\left(t\right)+{\eta }_H^{-}(t-T_H)]}^2-\frac{2}{\mu \kappa R_b}{\xi }_b^H\left(t\right)[T_b(\overline{\omega }+{\xi }_b\left(t\right))-T_b\left(\overline{\omega }\right)]\right\}}\hfill \end{array}$$

(98)

(see (74) and (91) for comparison). The novelty here is the pseudo-quadratic form

$$\begin{array}{c}{\displaystyle \mathcal{F}({\xi }_b^H,{\xi }_b)={\gamma }_0{c}_v^{\prime }{\left({\xi }_b^H\right)}^2-\frac{2{\gamma }_0}{\mu \kappa R_b}{\xi }_b^H[T_b(\overline{\omega }+{\xi }_b)-T_b\left(\overline{\omega }\right)]+{\xi }_b[T_b(\overline{\omega }+{\xi }_b)-T_b\left(\overline{\omega }\right)]+}\hfill \\ {\displaystyle +c_v{\xi }_b^2={\gamma }_0{c}_v^{\prime }{\left({\xi }_b^H\right)}^2-\frac{2{\gamma }_0}{\mu \kappa R_b}·\frac{T_b(\overline{\omega }+{\xi }_b)-T_b\left(\overline{\omega }\right)}{{\xi }_b}{\xi }_b{\xi }_b^H+}\hfill \\ {\displaystyle +\left[c_v+\frac{T_b(\overline{\omega }+{\xi }_b)-T_b\left(\overline{\omega }\right)}{{\xi }_b}\right]{\xi }_b^2.}\hfill \end{array}$$

(99)

Taking into account the sector condition (25) or (29), $\mathcal{F}({\xi }_b^H,{\xi }_b)>0$ provided the following inequality holds

$${\displaystyle {\gamma }_0{c}_v^{\prime }\left[c_v+\frac{T_b(\overline{\omega }+{\xi }_b)-T_b\left(\overline{\omega }\right)}{{\xi }_b}\right]>\frac{{\gamma }_0^2}{{\left(\mu \kappa R_b\right)}^2}{\left[c_v+\frac{T_b(\overline{\omega }+{\xi }_b)-T_b\left(\overline{\omega }\right)}{{\xi }_b}\right]}^2}$$

(100)

A rather elementary discussion, based on the 2nd degree trinomial properties shows that by choosing

$${\displaystyle {\gamma }_0>{\left(\mu \kappa R_b\right)}^2\frac{L+c_v}{L^2}{c}_v^{\prime },}$$

(101)

the strict positiveness of $\mathcal{F}({\xi }_b^H,{\xi }_b)$ is ensured provided (25) holds. Consequently $\mathrm{d}{\mathcal{V}}^{★}\left(t\right)/\mathrm{d}t\le 0$ and the stability of the composite system (70),(87) in the metrics induced by the Lyapunov functional defined by (97) follows.

For the asymptotic stability we shall turn again to the Barbashin-Krasovskii-LaSalle invariance principle. We examine firstly the asymptotic stability of the difference subsystem of the composite system: from (70) and (87) we can see that the two difference systems which define asymptotically stable difference operators do not interact. The characteristic equation of the difference operator is

$${\displaystyle \left(1+\frac{q-c_l}{q+c_l}{\mathrm{e}}^{-2\lambda T}\right)\left(1+\frac{q_H-c_l^{\prime }}{q_H+c_l^{\prime }}{\mathrm{e}}^{-2\lambda T+H}\right)=0}$$

(102)

and has its roots located on the vertical lines in ${\mathbb{C}}_{-}$

$${\displaystyle \Re e\left(\lambda \right)=\frac{1}{2T}ln\left|\frac{q-c_l}{q+c_l}\right|<0,\Re e\left(\lambda \right)=\frac{1}{2T_H}ln\left|\frac{q_H-c_l^{\prime }}{q_H+c_l^{\prime }}\right|<0.}$$

(103)

Next, the set where $\mathrm{d}{\mathcal{V}}^{★}\left(t\right)/\mathrm{d}t=0$ is defined in both (76) and (93). In this set - the kernel of the derivative function (98) - the equations of the composite system are defined by (77) and (94) which are decoupled. It follows from the previous development - done for each of these systems separately - that the only invariant set contained in the kernel of the derivative function (98) is the equilibrium at the origin. Asymptotic stability follows from the equilibrium at 0 for (70),(87). From the one-to-one correspondence between the solutions, the same property of the asymptotic stability for the equilibrium at the origin holds for the composite system of the Riemann invariants (65), (82) and for the basic composite system of the vibrations dynamics (62),(80). These results are summarized in the following

Theorem 7.

Consider the composite system of the coupled torsional and axial vibrations dynamics of the drillstring defined by (62),(80), together with the associate composite systems of the Riemann invariants (65),(82) and of the coupled delay differential and difference equations (70),(87) under the assumptions of Theorems 4 and 6. Then the equilibrium at the origin of each aforementioned system is globally asymptotically stable.

Obviously this theorem integrates the results of Theorems 4 and 6, strongly relying on Theorems 3 and 5.

6. Conclusions and perspective research

It was considered the dynamics of the coupled torsional and axial vibrations of the drillstring as obtained by making use of the Hamilton variational principle. A careful listing of the acting forces and torques allowed to complete the models in [36] such that the resulting difference operators of the associated systems of functional differential equations with deviated argument be asymptotically stable. As a consequence, after obtaining uniform stability via the energy like Lyapunov functional, it became possible to obtain asymptotic - even exponential - stability for the associated system of functional differential equations by applying the Barbashin-Krasovskii-LaSalle invariance principle. This property was then projected back on the initial system via the representation formulae for the Riemann invariants of the problem. The paper focused on the dynamics of the axial vibrations since this dynamics was forced by the coupling with the system of the torsional vibrations, the forcing signal approaching zero asymptotically as a consequence of the asymptotic stability of the torsional vibrations system. However, instead of using the results on perturbed dynamical system from e.g. [43,44,45,46,47,48,49,50], there was considered the coupled torsional/axial vibrations dynamics and a Lyapunov functional combined from those used separately. Once again asymptotic stability was obtained by applying the Barbashin-Krasovskii-LaSalle invariance principle for the overall system.

From our point of view the methodology applied in this paper is applicable to more complicated but possibly more adequate modeling containing e.g. nonlinear dependencies of the torques and forces on the state variables and other models of the drillstring viewed as a distributed parameter bar/beam/shaft [23,24]; the choice of the model depends on the choice of the elastic potential energy [25]. Moreover, the aforementioned aspects can be useful in the analysis of other mechanical systems.