Energy flow characteristics of friction-induced nonlinear vibrations in a water-lubricated bearing-shaft coupled system

Based on the energy flow theory of nonlinear dynamical system, the stabilities, bifurcations, possible periodical/chaotic motions of nonlinear water-lubricated bearing-shaft coupled systems are investigated in this paper. It is revealed that the energy flow characteristics around the equlibrium point of system behaving in the three types with different friction-paramters. (a) Energy flow matrix has two negative and one positive energy flow factors, constructing an attractive local zero-energy flow surface, in which free vibrations by initial disturbances show damped modulated oscillations with the system tending its equlibrium state, while forced vibrations by external forces show stable oscillations. (b) Energy flow matrix has one negative and two positive energy flow factors, spaning a divergence local zero-energy flow surface, so that the both free and forced vibrations are divergence oscillations with the system being unstable. (c) Energy flow matrix has a zero-energy flow factor and two opposite factors, which constructes a local zero-energy flow surface dividing the local phase space into stable, unstable and central subspace, and the simulation shows friction self-induced unstable vibrations for both free and forced cases. For a set of friction parameters, the system behaves a periodical oscillation, where the bearing motion tends zero and the shaft motion reaches a stable limit circle in phase space with the instant energy flow tending a constant and the time averaged one tending zero. Numerical simulations have not found any possible chaotic motions of the system. It is discovered that the damping matrices of cases (a), (b) and (c) respectively have positive, negative and zero diagonal elements, resulting in the different dynamic behavour of system, which gives a giderline to design the water-lubricated bearing unit with expected performance by adjusting the friction parameters for applications.


Introduction
Water-lubricated bearings have been more and more used in the marine and pump industries to eliminate the pollution of metal bearings lubricated by oil and grease as well as to increase efficiency and reliability of marine propulsion systems. The first effective water-lubricated bearing made of lignum vitae was invented by John Penn in the 1840s while man-made plastic bearings and natural rubber bearings appeared in the 1920s, which had gradually predominated in pumps, naval and many commercial ships by the 1940s. The need for improved water-lubricated bearings was greatly recognised in 1942 when a number of U.S. ships suffered extensive combat damage at the Battle of Midway, after which, USA navy department began to carry out an extensive research on water-lubricated bearings and obtained a series of achievements and the only military specification Executive Editor: Li-Feng Wang 1 3 of the water lubricated rubber bearings was created by MIL-DTL-17901C (SH) in 2005 [1]. Since the friction performance plays a key role in bearing vibration or noise, this specification especially describes the requirements for the test rig and scope of friction coefficients of bearing synthetic rubber facings.
Stern-tube bearing noise is associated with both the quantity of friction and the slope of friction-velocity curve [2], which involves bearing materials. To reduce bearing noises, many historical experimental and analytical investigations were reported [3][4][5]. For examples, Krauter [3,4] focused on a three-degree-of-freedom (3-DOF) system to incorporate with the essential elements of the friction and vibration phenomena in their experimental model, in which a linear analysis was considered to determine the instability conditions of the system that may lead to frictional vibrations. Sinou et al. [6] proposed a numerical model incorporating realistic laws of local friction in base of experimental results to characterize the dynamics of a lubricated system and to study its complex global responses triggered by the local interfacial behaviour. Graf and Ostermeyer [7] have shown how the stability of an oscillator sliding on a belt will change, if a dynamic friction law with inner variable is considered instead of a velocity-dependent coefficient of friction, which demonstrates the unstable vibration can even be found in the case of a positive velocity-dependency of friction coefficient. Baramsky et al. [8] contributed the analytical and experimental results for the occurrence of friction-induced vibrations during tightening of bolted joints of the most used machine elements. Niknam and Farhang [9] employed a two-degree-of-freedom single mass-on-belt model to study friction-induced instability due to mode-coupling, in which numerical analyses are used to tackle the effect of three parameters related to belt velocity, friction coefficient, and normal load on the mass response. Ghorbel et al. [10] proposed a minimal 2-DOF disk brake model to investigate the effects of different parameters on mode-coupling instability, which considers self-excited vibration, gyroscopic effect, friction-induced damping, and brake pad geometry.
In designs of this type of berings, two research directions have been developing: one aims to obtain good friction characteristics using suitable geometric parameters with different facing layers and lubrications [11][12][13][14], and another is to create new bearing materials [15][16][17][18]. Recently, Smith [19] presented a paper on the design and lubrication of waterlubricated, rubber, cutlass bearings systems, which gave a methodology to predict the minimum film thickness between the journal and heavily loaded stave, including a 3D finite element (FE) model to predict generated pressures in cutlass bearings, and comparing with experimental data.
From the point view of sciences, the water-lubricated bearings system is a typical nonlinear system concerning friction-induced vibration noises, so that for good noise performance of this type of bearing systems, it is essential to investigate its integrated nonlinear characteristics. In the area on friction induced vibrations, Ibrahim [20,21] contributed two important review papers, of which the first one concerns the mechanics of contact and friction and the second one discusses the dynamics and modeling of friction induced vibrations. Both of papers provide the comprehensive account of the main theorems and mechanisms developed in the historical literatures concerning friction-induced noise and vibrations with a quite wide references listed. Furthermore, Ding and Zhai [22] presented a review paper on the advances of friction dynamics in mechanics system, in which the common used friction models, friction related self-excited vibrations and their controls were presented. The interested readers may refer these three review papers for more information and references in this area.
More directly linking the topic of this paper, Simpson and Ibrahim [2] investigated the dynamic behaviour of nonlinear friction-induced vibration in water-lubricated bearings based on the traditional method dealing with nonlinear dynamical systems. Their paper developed an analytical nonlinear, two-degree-of-freedom model, which emulates the stern-tube bearing. Typical friction-speed curves were adopted based on the experimental results of Krauter [23]. The stability conditions of equilibrium were predicated. In the unstable cases, the nonlinear response behavior was examined by using numerical integration of the coupled equations of motion. The dependence of the relative sliding speed on time and its effect on the friction force was included in the numerical simulation. A control parameter that combines the influence of the normal load and the slope of the friction-speed curve was used to construct a bifurcation diagram which separated different response characteristics, including squeal, limit cycle, and stable regions.
To reveal nonlinear behavior of nonlinear dynamical systems, such as as discussed in the books by Guckenheimer and Holmes [24], Thompson and Stewart [25], Chen et al. [26] and Liu and Chen [27], Xing [28] developed a genralised energy flow theory to investigate nonlinear dynamical systems governed by ordinary differential equation in phase space. Important nonlinear phenomena such as, stabilities, periodical orbits, bifurcations and chaos can be revealed from the energy flow behavors of the systems [29]. The aim and motivation of this paper are to investigate the 2-DOF nonlinear model of the stern-tube bearing proposed in Ref. [2] to reveal its energy flow characteristics governing its nonlinear dynamical behavours, which is also as an example to further develop energy flow approach based on two scalar varaibles, genralised potentional and kinetic energies, to tackle complex nonlinear dynamical systems, especially in multi-dimensional phase spaces.

Dynamic modelling and equations of water-lubricated bearings
The marine propeller shafts are supported in stern tube bearings, which are lubricated by water and are almost all grooved, of which the scketched structure is shown in Fig. 1. The 2-DOF analytical model used in Ref. [2] is shown in Fig. 2, where the mass, stiffness, and damping coefficient of the bearing (subscript 1) and the flexible shaft (subscript 2) are denoted by m 1 , m 2 , k 1 , k 2 , c 1 and c 2 , respectively. Here, the forces F 1 and F 2 denote possible prescribed external forces applied to the bearing and the shaft, respectively, while the force F f denotes the friction force acted to the bearing by the shaft, which is assumed in the viscous type propotional to the relative sliding speed of the friction pair with a viscous damping coefficient b , so that F f = b(̇x 2 −̇x 1 ).

Viscous damping coefficient
Historically, the friction force was formulated by the friction model developed by Smith et al. [5] and used by Simpson and Ibrahim [2] in their 2-DOF analytical model, that is where N denotes the normal bearing load; a is a constant with units of inverse velocity; 0 and 1 are the static and dynamic coefficients of friction, respectively; R represents the disk radius and ̇ is the torsional velocity of the disk relative to the bearing. Considering only oscillatory motions of both the shaft and bearing, and representing the motion by the linear displacement, we can denote the relative sliding speed of the friction pair in the form Substituting Eq. (2) into Eq. (1) and then taking the derivative of the resultant equation with respect to the relative sliding speed, we obtain the visous damping coefficient as which, when the function e −av is expanded into a Taylor series, becomes

Dynamic equations of the friction pair system
Using the model shown in Fig. 2 and the second Newton's law, we obtain the following two equations of motion of the system which can be rewritten in a matrix form where the mass, damping and stiffness matrices as well as the displacement and force vectors are resepctively defined as follows Introducing the nondimensional parameters

Bearing housing
Rubber stave Gap filled with water Rotor shaft  where i denote dimensionless damping coefficients, 21 ratio of circular frequencies, i natural frequencies of the bearing and the shaft subsystems, 12 mass ratio, y i dimensionless displacements, dimensionless time, dimensionless friction coefficient, −1 dimensionless normal bearing load and the dimensionless value of constant . Using these dimensionless parameters, the dynamic equation of system can be written as where For the approximation as shown in Eq. (4) adopted by Simpson and Ibrahim [2], the parameter b and its derivatives are given by the following approximated formulations implying that the parameter b is approximated by a quadratic function of the relative velocity (z 2 − z 1 ) , its first derivatives are linear, and the second derivatives are constants.

Equations in phase space
Using the state variables y and z , we can rewrite Eq. (9) in the form of the phase space, i.e. of which, the vector (9) mÿ + ċy + ky =f , , c =ĉ +c,ĉ = 2 1 0 0 is called the generalised force of system. Here, J denotes the coefficient matrix of the state variable vector [y T z T ] T , which involves the nonlinear damping parameter b . The solution of Eq. (11) with a prescribed initial condition gives an orbit starting from the position r 0 in the phase space as shown in Fig. 3. The parameters of system, such as k 1 , c 1 and c 2 , can be functions of normal load and temperature. Krauter [4] developed empirical relationships for these parameters Typical values of these parameters were selected for constant values of normal load and temperature, i.e. N = 120 N (27 lb) and T =35 ℃, by Simpson and Ibrahim [2] to investigate the nonlinear dynamical behaviour of system. We also use these chosen parameters to reveal the energy flow characteristics of system, and to rewrite Eqs. (11) and (12) in the form from which in association with Eq. (10) we obtain (13)

Equilibrium point and Jacobian matrix
The point at which F = 0 is called the equilibrium point of system. Obviously, if there is no external force ( f = 0 ), the origin (y = 0 = z) of phase space is an equilibrium point of system. By using Eqs. (15) and (16), the Jacobian matrix of the system can be derived as Here the gradient operator ∇ is defined by Substituting Eq. (15) into Eq. (17), we finally obtain the Jacobian matrix Here the matrices J 0 and J b are defined by Eq. (15).

Generalised potential energy and kinetic energy
We define two scalar variables: the generalised potential energy and generalised kinetical energy of system in the following forms. Generalised potential energy and its time average Generalised kinetical energy and its time average Physically, these two scalars may not practical potential and kinetical energies, therefore we use the word generalised to distinguish them with practical physical quantities. Geometrically, as shown in Fig. 3, the generalised potential energy equals the half square of the distance of a point to the origin of phase space, i.e.
while the generalised kinetical energy equals the half square of the tangent vector, the velocity, at a point, i.e.
Therefore, these two non-negative real scalars can be used to describe the postion and velocity at a point on the solution orbit of system in phase space, and to reveal related nonlinear dynamical behavours of system.

Energy flow equation and its time average
Pre-multiplicating Eqs. (13) and (15) by a row vector r T = y T z T and its initial value r T (0) , respectively, we obtain the energy flow equation of which, P denotes the power done by the generalised force of system, p is the power of the external forces and Ė represents the time change rate of generalised potential energy, called as the energy flow of system. Physically, the energy flow equation in Eq. (24a) implies that the energy flow of system equals the power done by the generalised force.
The energy flow equation in Eq. (24a) can be rewritten in the form Here Ẽ is a real symmetrical matrix, while Ũ is a real antisymmetrical matrix.

Zero-energy flow surfaces and equilibrium points
Generally, the energy flow of system is a function of time and the position of a point in phase space, which generates a scalar field called as the energy flow field of the nonlinear dynamical system. Equation defines a generalised surface or subspace in phase space, which is called as a zero-energy flow surface on which the energy flow vanishes. If an orbit of the nonlinear dynamical system is on a zero-energy flow surface, the distance of a point on the orbit is not changed with time.
Based on the theory of extreme values of a function and the geometrical meaning of generalised potential energyE , we can conclude that the characteristics of orbit at a point on the zero-energy flow surface as follows.
• Ë > 0 , the generalised potential energy takes a local minimum-extreme value at this point compared with the points in the local domain around it. Therefore, with time going, the value of E increases and the orbit backwards the origin of phase space. • Ë = 0 , it cannot determine that generalised potential energy takes a local minimum-or maximum-extreme value at this point, and higher time derivatives are needed to give a solution. • Ë < 0 , the generalised potential energy takes a local maximum-extreme value at this point compared with the points in the local domain around it. Therefore, with time going, the value of E decreases and the orbit towards the origin of phase space.
If no external force p = 0, there are the following three cases satisfying Ė = 0 in Eq. (26).
• Case 1:r = 0 , that represents the origin of phase space, at which the generalised potential energy is defined as zero; • Case 2:F = 0 , implying an equilibrium point of system, so that equilibrium points are on the zero-energy flow surface; • Case 3: P = 0, r ≠ 0 ≠ F, correseponding a generalised zero-energy flow surface.
Assume that r denotes a point on an orbit on a zeroenergy flow surface, P(t, r) = r T F = r i F i = 0 , and is an small orbit variation around r , generally, the variation of energy flow caused by the orbit variation can be approximated to the quantities of | | 2 in the form where the energy flow gradient vector p , the energy flow matrices E and E 1 , and spin matric U are given by Here, we have used Eqs. (10a) and (15) and noticed that the components F 1 and F 2 of the vector F are linear functions of y and z , so that their second derivatives vanish.
For an equilibrium point satisfying F(r) = 0 , Eq. (27a) becomes which is further reduced to for an equilibrium point r = 0, and so that E 1 = 0. At the equilibrium point r = 0 the second time derivative of the generalised energy in Eq. (24a) vanishes, so that we need to consider the variation of Ė around the point r = 0 to investigate the orbit behaviour. From Eq. (28a) and the geometrical meaning in Eq. (22), we consider the variation of energy flow Ė to determine the characteristics of orbit as follows.
implies the flow towards the zero-energy flow surface, while ΔĖ < 0 indicates the flow backwards the zeroenergy flow surface; implies the flow towards the sero-energy flow surface, while ΔĖ > 0 indicates the flow backwards the zeroenergy flow surface; • If the flows from both sides of the zero-energy flow surface tend to it, this surface is an attracting surface. • The local ability of equilibrium point r = 0 can be determined by the behaviour of energy flow matrix E as follows: Figure 4 shows a case where the orbit intersects at a point r on the energy flow surface. Since ΔĖ > 0 above the surface and ΔĖ < 0 under the surface, so that the flow along the orbit backwards this point and this point is an unstable point.

Energy flow matrix and energy flow characteristic factor
Both of matrices Ẽ and E are real symmetrical matrices, called as the energy flow matrices, the former is for the total energy flow in Eq. (24b), while the later for the incremental energy flow relative to the zero-energy flow surface in Eq. (27b). The real symmetrical matrix has its real eigenvalues I and corresponding eigenvector I satisfying the orthogonal relationships, such as for matrix Ẽ in Eq. (24b), we have These eigenvectors span an energy flow space in which the vector r can be represented by which, when substituted into Eq. (24b), gives Therefore, ̃I represents the energy flow variation caused by a unit disturbance ̃2 I = 1 in the I-th principal direction ̃I of the energy flow matrix. Therefore, we call definitely -negative, asymptotic stable, semi -definitely -negative, stable, definitly or semi -definitly -positive, unstable.
as the energy flow characteristic factors and we can conclude for a point of the orbit of nonlinear system.
• The positive, zero and negative value of the factors ̃I respectively implies the energy flow increase, unchanged and decrease caused by the disturbance ̃2 I in the I-th principal direction, from which the behaviour of energy flow increments at this point can be identified to judge the local dynamical behaviour, such as for local stability of orbits around an equilibrium point r = 0, we have • If its energy flow characteristic factors are not all seminegative or not all semi-positive, there will exist a small subdomain around this point in the phase space determined by Ė = 0 , from which a zero-energy flow surface can be obtained.

Spin matrices and periodical orbit
The matrices Ũ in Eq. (24b) and U in Eq. (27b) are two real skew-symmetrical matrices, called as spin matrices. For periodical orbits, it is neccesary that there exsist a time period T and its corrsponding closed orbit Γ such that the following integrations along the closed orbit hold, i.e.
where Ê and K are two positive constants. The stability of this periodical orbit can be determined by the values of the energy flow ΔĖ around each point of the orbit as discussed for the zero-energy flow surface. The time averaged kinetic energy can be re-written as Physically, this integration denotes the circulation integral of velocity field ̇r along the closed orbit and involves the skew-symmetrical spin matrix. To clarify this, we consider a 3-D vector field ̇y = f , the curlf , or ∇ × f at a point O is explained in Fig. 5. Here is a unit vector, the projection of the curlf onto is defined as the limit of a closed line integral along the curve C in a plane orthogonal to , that is negative , stable subspace, zero, central flow subspace, positive, unstable subspace. where A is the area enclosed by curve C. The curlf can be denoted in a tensor form [30,31] where e ijk is the permutation tensor. The vector curlf is a dual vector of a skew-symmetrical matrix U , spin matrix, satisfying the following relationship Therefore, a positive time-averaged kinetic energy implies there must be a non-zero spin matrix for periodical orbits.

Bifurcation and chaos
The energy flow of a nonlinear dynamical system is affected by the bifurcation parameters to reveal the birfurcation characteristics of orbits. For example, different parameters result in different equilibrium points, zero-energy flow surfaces and energy flow characteristic factors etc.
Xing [28,29] has discovered that a chaotic motion of nonlinear dynamical system can be considered as a periodical motion with an infinte long time period and the following integrations hold (34c) e ijk e krs f s,r This implies that the time averaged mechanical energy ⟨E + K⟩ tends a constant when the average time tends to infinite. Also, for a chaotic motion, flows are restricted in a finite volume, so that the space averaged rate of volume strain of phase space must not be positive, i.e.

Investigations of nonlinear
water-lubricated bearing system

Energy flow characteristic factors and vectors
For the approximated energy flow matrix E , using Eqs.

Bifurcation of zero-energy flow surfaces around the fixed point
The dimensionless parameters studied in Ref. [2] are listed as follows based on which the obtained energy flow characteristic factors around the fixed point are given in Tables 1, 2 Fig. 6a, and the computer simulating surface is drawn in Fig. 6b, while the zero-energy surface based on Eq. (37) for nonlinear case is given in Fig. 6c. In Fig. 6a, the three characteristic vectors ⃗ 2 , ⃗ 3 and ⃗ 4 respectively with the three factors span a subspace in phase space o − y 2 z 1 z 2 . This local linearized zero-energy flow surface is symmetrical to the coordinate plane o − ⃗ 2 ⃗ 3 and axis-symmetrical to the axis o − ⃗ 4 , so that a plane perpendicular to axis o − ⃗ 4 intersects the zero-energy flow surface in an elliptical curve of semi-axis | |̃4 | | R 3 > | |̃4 | | R 2 , as shown in Fig. 6 for the ones at | |̃4 | | = 1 . The plane o − ⃗ 3 ⃗ 4 intersects the zero-energy flow surface in the straight lines of ̃3 = ±̃4R 3 . Based on Eq. (42a), at a disturbed point outside the zero-energy flow surface on the plane ̃4=const., the energy flow Ė < 0 , so that it intends to reduce the distance to the origin of space, implying the flow towards the zeroenergy surface. For a disturbed point inside the zero-energy flow surface on the plane ̃4=const., the energy flow Ė > 0 and it intends to increase the distance to the origin of space implying the flow also towards the zero-energy surface. Therefore, the zero-energy surface is an attracting surface. In this case, from Eq. (10a), the damping matrix of system at the fixed point takes with the positive diagonal damping coefficients in the system, so that the motions caused by initial conditions will be damped.
By using the MATLAB program in Ref. [28] and the nonlinear friction parameter in Eq. (37b), the results of numerical simulations with the initial conditions y = 10 −6 1 1 1 1 T for this case are given in Figs. 7, 8, 9. Figure 7 gives the time histories of displacement and velocity of the bearing and the shaft, which shows the damped oscillations. Figure 8 presents the phase diagram y(1)-y (3) of bearing unit and the one y(2)-y(4) of shaft unit, as well as the time histories of generalized potential energy and the distance of phase point to the origin, which indicates the orbit points in both phase diagrams move from outside to the origin, and the generalized potential energy and the distance of orbit point to the origin tend to zero with time increasing. Figure 9 shows the instant energy flow behaving a damped oscillation and the time history of phase space volume strain decreasing with time. ,

3
The sketched zero-energy flow surface and energy flow variations around the equilibrium point is shown in Fig. 10a, while the computer drawn surface based on Eq. (44a) is given in Fig. 10b. The one for the nonlinear equation is like to the one shown in Fig. 6c, so neglected. The local linearised surface is symmetrical about the coordinate plane o − ⃗ 3 ⃗ 4 and axis-symmetrical to the axis o − ⃗ 2 , so that a plane perpendicular to axis o − ⃗ 2 intersects the surface in an elliptical curve of semiaxis | |̃2 | | R 3 > | |̃2 | | R 4 , as shown in Fig. 10a for the ones at | |̃2 | | = 1 . The plane o − ⃗ 2 ⃗ 3 intersects the surface in the two straight lines of ̃3 = ±̃2R 3 . At a disturbed point outside the zero-energy flow surface on the plane perpendicular to the axis o − ⃗ 2 , the energy flow Ė > 0 by Eq. (42a), so that this disturbance intends to increase the distance to the origin of space, implying the flow backwards the zero-energy surface. For a disturbed point inside the zeroenergy flow surface, the energy flow Ė < 0 and the disturbance intends to reduce the distance to the origin of space implying the flow also backwards the zero-energy surface. Therefore, this zero-energy surface is a divergence surface. In fact, for this case, the damping matrix of the system in Eq. (10a) becomes that shows the negative diagonal damping coefficients for both bearing and shaft, so that the motions caused by initial conditions will be divergence.
The results of numerical simulations with the initial conditions y = 10 −6 1 1 1 1 T for this case are given in Figs. 11-13. Figure 11 gives the time histories of displacement and velocity of the bearing and shaft, which shows the   Figure 12 presents the phase diagram y(1)-y(3) of bearing and the one y(2)-y(4) of shaft, as well as the time histories of generalized potential energy and the distance of phase point to the origin, which indicates the orbit points in both phase diagrams move far away to the origin, and the generalized potential energy as well as the distance of orbit point to the origin increase. Figure 13 shows the instant energy flow and the time history of phase space volume strain, of which both positively increase very fast at about non-dimension time 350. In this case, since there exsists a zero-energy flow characteristic factor, the dynamic behavour of the system may not be determined only based on linearised approximation as discussed below. In this case, the central energy flow investigation by considering higher-order terms is needed [28]. The zero energy flow serface around the origin becomes corresponding two orthogonal planes perpendicular to the coordinate place o − ⃗ 2 ⃗ 4 with the intersect lines and ̃2 = −̃4, respectively, as shown in Fig. 14. In the domain containing axis o − ⃗ 2 , the energy flow Ė < 0 , so that the orbit disturbance along o − ⃗ 2 axis will be towards the zero-energy flow surface, while in the domain containing o − ⃗ 4 axis, Ė > 0 , the disturbance along o − ⃗ 4 axis will be backwards the zero-energy surface. Therefore, the dynamic behaviour of the system cannot be predicted only based on linear approximation. The theoretical analysis relies on central energy flow theorem, which is more complex, so that (45a) 0.512142̃2 2 = 0.512142̃2 4 ,̃2 = ±̃4, here we will reveal the solution by numerically simulating the case with above parameters. Also, from linearised approximation, the summation of diagonal elements of the energy flow matrix E vanish, so that its time change rate of phase volume strain approximately vanishes, i.e.

Case
implying the volume of phase diagram should not change, which might not be true for the nonlinear equation with no linearazation. Now, we show the numerical simulation results to explain above discussion. The damping matrix at fixed point takes value of which the diagonal damping coefficients of system vanish. Figure 15 shows the time histories of displacement and velocity for both bearing and shaft units, from which it has been found that before the time 9.4 × 10 4 s, the force on shaft f = 0 0 0f T • Case (a) when the energy flow matrix at the equilibrium point has two negative and one positive characteristic factor, the local zero-energy flow surface in the small domain around the equilibrium point behaves an attractive surface, and the system is stable to initial disturbance. We also simulated some cases reported in Ref. [2], which further confirm this conclusion is correct. For example, the three cases with = 2.0, 1 = 0.1116, 0.1111, 0.111245, studied based on the friction force approximation in Eq. (37c) in Ref. [2] were reported the stable case for 1 = 0.1116 but unstable cases for the last two values of 1 = 0.1111, 0.111245 . We re-check these three cases by using the friction force Eq. (37b), and the simulation  Tables 1, 2 and 3, respectively, and the equilibrium point is stable. • Case (b) when the energy flow matrix at the equilibrium point has one negative and two positive characteristic factors, the local zero-energy flow surface in the small domain around the equilibrium point behaves a divergence surface, and the system is unstable. • Case (c) when the energy flow matrix at the equilibrium point has one zero factor and two same absolute value factors, one negative and another positive, the local zeroenergy flow surface in the small domain around the equilibrium point behaves divergence in a sub-domain and attractive in another sub-domain as shown in Fig. 14. The dynamical behaviour about this case requires to undergo a higher order approximation analysis. The numerical simulation result shows a divergence self-excited oscillation, which might be miss-judged as a stable oscillation if the simulation time is not enough large. For example, Fig. 18, cut out from Fig. 15 ( = 2.0, 1 = 0.111 ), shows the enlarged time histories of displacement and velocity of both bearing and shaft in the initial time before the simulation breaking, which seems a stable picture but unstable indicated by Fig. 15.

Periodical oscillation (f = 0)
From Sect. 2.6, we know that for a possible periodical orbit, a necessary condition is that the spin matrix U of the system must not vanish. From Eqs. (15), (19) and (27b), we obtain the spin matrix at the equilibrium point  Figure 19 shows that the displacement and velocity of bearing decrease with time, while the shaft ones are in a periodical motion. The phase diagram of bearing tends to the origin while the one of shaft shows a closed orbit, and the generalised potential energy as well as the distance c = 0.0012 0.2220 0.2220 0 .
of phase point to the origin behave modulated oscillation without obviously amplitude variations, as shown in Fig. 20. The instant energy flow curve behaviours a stable oscillation in Fig. 21a, and its time averaged one tends to zero with average time increasing as shown in Fig. 22. The phase space volume train descreases with time shown in Fig. 21b. According to the energy flow theory discussed in Sect. 2.6, this is a stable periodical motion of the system excited by initial conditions. We have simulated the cases with parame t e r s ( = 2.0 1 = 0.1111, 2 = 0.008325) a n d ( = 2.12 1 = 0.119725, 2 = 0.0009) [2], of which the energy flow matrix has two negative and one positive energy flow characteristic factors and both systems also behave the similar stable periodical oscillations as the one presented herein.

Forced oscillations (f ≠ 0)
We respectively add an excitation force of amplitude 10 −6 , frequency 0.5 and phase angle 0, i.e. f = 10 −6 cos(0.5 ), on the bearing mass and the shaft mass to investigate the forced dynamic response of the system. The simulated results are given as follows.  Fig. 23 tend stable forced modulated oscillations, the generalised potential energies and the distances of phase points to the origin tend stable pictures and the phase diagrams tend periodical orbits in Fig. 24, as well as the instant energy flows and the phase space volume strains in Fig. 25 also tend stable curves. Figure 26 confirms that the time-averaged energy flows tend zero with average time increasing. These typical results demonstrate that the forced vibrations are periodical motions of the system although the different curve styles due to different force added potions. Since the positive damping of the system, the free vibration components due to initial conditions are gradually damped, and a stable forced vibration is obtained. However, since the system is nonlinear, the frequency of the forced vibration does not equal the force frequency, so that the response curves behave modulated pictures compositing of two frequencies response. The forced vibration in this case is like the free vibration discussed for case 3 in Sect. 3.3, as shown in Fig. 29, the dynamic responses of displacements and velocities of both bearing and shaft increase until about time 3.3 × 10 4 when the responses sharply increased and simulations stopped.
To show the response curve details, Fig. 30 provides local enlarged curves before breaking time in Fig. 29, from which it is observed the dynamic responses behave the modulated vibrations before the simulation breaking. The curves of generalised potential energy, the phase diagram, etc. also show the characteristics, so they are neglected herein.

Chaotic motions
With different friction parameters, we have tested several cases for free vibrations induced by given initial disturbances and for forced vibrations excited by external forces acted on the bearing and shaft masses in order to see if chaotic motions may be found. The simulations results have not shown any chaotic motions of this system.

Conclusion
Following the detailed describption of this investigation given above, we may conclude the main contributions as follows. The paper is a first one to study a nonlinear waterlubricated bearing-shaft coupled system using the energy flow approach, which practices and further confirmed an efective energy-flow means to tackle varous complex nonlinear systems and to investigate their nonlinear characteristics: stability, periodical motion, birfurcation and possible chaos.The results obtained for the investigated system have revealed that the dynamic behaviour about its equilibrium point is fully determined by the energy flow matrix and its characteristic factors. Three cases with different energy flow charactristic factors affected by the friction-parameters are discovered: (a) two negative and one positive factors, constructing an attractive local zero-energy flow surface, in which free vibrations show damped modulated oscillations allowing the system returning to its equlibrium state, while forced vibrations show stable oscillations; (b) one negative and two positive factors, spanning a divergence local zeroenergy flow surface, so that the both free / forced vibrations are in divergence resulting an unstable system; (c) one zero and two opposite factors, constructing a local zero-energy flow surface dividing the phase space into stable, unstable and central subspace, and the simulation shows friction selfinduced unstable vibrations in both free and forced cases. For a set of friction parameters, the system behaves a periodical oscillation, in which the bearing motion tends zero and the shaft motion tends a stable limit circle in phase space with the instant energy flow tends a constant and its time averaged one tends zero. Numerical simulations have not found any possible chaotic motions of the system. It is discovered that the damping matrices of cases (a), (b) and (c) respectively have positive, negative and zero diagonal elements, dominating the dynamic behavour in different cases of the system, which provides an approach to design the water-lubricated bearing unit with expecting performance in marine engineering applications by choosing suitable friction parameters to generate the expected damping matrix.