The determining torsional vibration damping coefficients algorithm for computing marine shafting’s vibrations

411
The determining torsional vibration damping coefficients algorithm for
computing marine shafting’s vibrations
Ass. Prof. Nguyen Manh Thuong
Vietnam Maritime University, Email: thuongmt29@gmail.com
Abstract The most difficulty in the computations of ship shafting torsional vibrations is to determine
such torsional vibration damping coefficients as those in main engines and shafting structures, or those
of hydraulic resistances on propellers. The difficulties are due to the complex nature of the hydro
dynamical process, that it is difficult to come up with an analytical method to determine these quantities.
Even direct measurements of such quantities are difficult to make in the practice.
Currently in Vietnam, these data usually are given by the manufacturers that doesn’t make sure about
the accuracy, or they are calculated based on the semi-empirical formulae those are old and may be not
appropriate for the modern types of engines and equipment. Therefore, the aim of the author is to
determine these quantities from the results of torsional vibration measurements during the new built ship
tests. Based on the analysis of results obtained for shafting’s types with different features (for example,
characteristic types of engines, propellers, shaft materials ...) and combined with the theoretical analysis
we would generalize calculations of the above quantities. It allows to obtain the new calculating
formulae for these quantities more certain and convenient for basic design calculation of ship shafting.
The paper presents a theoretical basis and algorithms of determining the torsional vibration damping
coefficients of the engine structure, shaft material and hydrodynamic damping coefficient to propellers
from results of vibration measurements. By programming with the help of Symbolic Math Toolbox
package in Matlab software we will establish an analytical equation for amplitudes of elastic torque
( )M k caused by torsional vibrations in an arbitrary shaft section (where perform measurements
conveniently). In this equation it contains m unknown variables, supposed to be unknown damping
coefficients for each given shafting with given size and mass parameters and a finite number of masses.
Therefore, to identify these m unknown variables, we must have at least m values of torque amplitude
corresponding to m values of k measured in a position of shafting to make a system of m equations.
However, the system of equations is not linear so to facilitate the calculations we need to have more
than m values of the measured amplitudes. The program would indicate how many values of the
measured amplitudes should have, and after having the measurement results, it will determine the
unknown damping coefficients.
To illustrate, the article presents the example calculation for ships series B 170-V designed by Polish
Architecture Institute and built in Ha Long shipyard. The data on the shafting and the results of vibration
measurements obtained from the documents submitted by Ha Long shipyard and the Designer to GL
registry for approval. The obtained results of calculation will explain the reason for deviation between
real vibrations and computed vibrations or show the appropriateness of assumptions about the damping
coefficients used in computations of torsional vibrations. Based on this profile, we can be able to
improve the vibration calculation model.
Keywords: torsional vibrations, damping coefficient, relative and absolute damping coefficients, elastic
torque.
1. Introduction
As we know, the kinds of resistances acting on the shafting’s vibrations in general are not constant, but
non-linearly depend on amplitudes and velocities of vibrations. For example, the resistance in the shaft’s
material (relative resistance) or hydraulic resistances on propellers (absolute resistance) is the
exponential functions of torsional deformation amplitudes or torsional vibration amplitudes of propellers
(for example, the resistant work to propellers per a vibration cycle is proportional to square of the
412
vibration amplitude [1]). Therefore here we should solve systems of nonlinear differential equations. It
is the reason why such complicated resistant models are used to calculate only resonant vibrations,
assuming that forms of forced vibrations are similar to forms of free vibrations [1, 3].
But currently, many Design Agencies and Registries have accept assumptions, supposing these damping
coefficients are constants and resistant torques proportional to the speed of deformations. The damping
coefficients of the water on propellers accord to the Archer formula [2].
Figures 1a and 1b show calculated by the Polish design agency and measured torques of elastic
deformations at intermediate shaft of B-170-V ship and figure 2 shows the calculation results of the
authors of this paper. As seen on the pictures, the frequency curves of calculations and of actual
measurement are rather identical. Differences in quantities are interpreted as errors caused by using
simplifying assumptions about resistances.
Figure 2 Stress amplitude/ rotation speed at intermediate shaft calculated by the
author
In spite of the fact that there are differences between computations and measurements, we may see that,
computation results with the above simplifying approximate assumptions meet the purposes of torsional
vibration calculation for basic design stage and the requirements of the Registry for torsional vibration
calculations. But essential problem posed now is how to determine these damping coefficients. How do
Figure 1a Stress amplitude/ rotation speed at
intermediate shaft calculated by the Polish
designers [2]
Figure 1b Stress amplitude/ rotation
speed at intermediate shaft measured
[2]
413
we design without data given by manufacturers? What is to be based on to confirm results of design
vibration calculations be reliable?
We may infer that, the determination of the above damping coefficients still is the problem to be solved
even for the world developed shipbuilding industry because the foreign registry offices still require
testing torsional vibration measurements for new designed ships.
This paper will present algorithms determining damping coefficients having had the torsional
measurement results when testing new build for the use in the future to study these types of resistances.
Indeed, after determining the damping coefficients for many different ship shaftings, we will have more
experience and understanding of the factors causing resistances that may help us to establish semi-
empirical formulas to determine the coefficients of resistances in engines and in the shaft material or
resistances to propellers by water ... to serve torsional vibration calculations in design stage.
2. Establish algorithms determining the damping coefficients
Let’s consider an actual shafting as a torsional vibration system of n discrete masses, linking together
by elastic elements without weights as shown in Figure 3 [1, 2]. The vibration driving forces are forces
from cylinders acting on engine crank shafts and of water on propeller. The main forms of resistance
included are those in the structure of piston-cylinder/connecting rod/crank shaft, and those in shaft
material and of water on propellers. Designating angle of torsion deformation of i-th mass as φi, forcing
torque on that mass as Mi, inertial mass moment of the i-th mass as Ii, absolute damping coefficient- ai,
relative damping coefficient between two masses i and i+1 as bi,i+1 and stiffness between them- Ki,i+1,
we have dynamical equation of i_th mass as [3]:
1, 1 1 1 1, , 1 1, , 1 , 1 1 , 1 , 1( ) ( )i i i i i i i i i i i i i i i i i i i i i i i i i ib K I a b b K K b K M    
where: i = 1,,n. (1)
Figure 3 Torsional vibration model of B_170 V, main engine - 6 RTA 62 U
For stable vibrations, torsional angle of i-th mass is equal sum of harmonic vibrations:
, (2)
and the exciting moments, being cyclic, can be presented by Fourier’s series:
. (3)
Replacing (2) and (3) into (1) we will get an equation system for each harmonic order k, after equating
coefficient of sin and cos we will get a system of 2×n algebraic equations as [3]:
, (4)
    , ,
1
( sin cos )
m
i i k i kS k t C k t
     , , ,
1 1
sin( ) ( sin cos )
m m
s c
i i k i i k i kM M k t M k t M k t
 k k kF A M
414
where:
- column vector of exciting torque amplitudes at
harmonic order k;
- column vector of torsion angle amplitudes according to
harmonic order k;
Fk = [fi,j ] – square matrix 2×n, its elements fi,j are functions of kω, stiffness of shafts Ki,i+1, inertial
moment Ii and of damping coefficients ai, bi,i+1 [4].
In calculations forced vibrations, the elements fi,j of Fk are determined if damping coefficients are given
and torsional vibration amplitudes of masses due to harmonic k are computed by:
. (5)
Then elastic torsion torque amplitude at order k in shaft between masses i-th and (i+1)_th is calculated
by:
. (6)
If supposing we measure the elastic torsional torque in a shaft Ei,i+1 (see Figure 4), then according to (6)
we can determine the amplitude of the relative torsion deformation
, 1i i between two masses i and
i+1:
. (7)
Figure 4 Stress/time at the intermediate shaft at different rotation speeds [2]
On the other side
, 1i i can be calculated by equation (5) as following. Having:
  1,1...2 ,1...2 1, , .det( )i n i n k i k i k kM S S F   
  1,1...2 ,1...2 1, , .det( )n i n n i n k i k i k kM C C F   
Therefore:
 2 22 2 2, 1 1,1...2 ,1...2 1,1...2 ,1...2det ( ) .i i k i n i n i n i n kF M     ,
(8)
1, 2, , 1, 2, ,=[ ,. ..., , , ...., ]
s s s c c c T
k k k n k k k n kM M M M M M M
 1, 2, , 1, 2, ,[ , ,..., , , ,..., ]
T
k k k n k k k n kA S S S C C C
 1k k kA F M
2 2
, 1 , 1 , 1, , 1, , 1 , 1( ) ( )
k
i i i i i k i k i k i k i i i iE K S S C C K
 , 1 , 1 , 1/
k
i i i i i iE K
415
where: ; ; 1,1...2n i n  ; ,1...2n i n  are lines i, i+1 and n+i+1, n+i of matrix
with elements determined as: - where - algebraic complement of the element (j,i) of Fk.
, 1i i can be calculated according to (7) if the elastic deformation torque is measured.
Thus, suppose there are m values of damping coefficients aj and bj, j+1 are unknown and should be
defined. Then equation (8) will contain m unknown variables (they are damping coefficients in Fk). In
principle, in order to determine the m these variables, we should get m values of elastic torsional torque
in shaft between masses i, i+1 (see F. 4) in order to make m algebraic equations to determine m unknown
values of aj and bj,j+1. However, the system of equations obtained is the type of multi-hidden variable
and high order, so that to be able to solve it, we should make more than m equations to obtain a system
of linear equations with the new hidden variables, those are combination of productions of unknown as::
To illustrate, here we consider the ships of series B 170 V. Its torsional vibration model is as shown in
F.3 having the following characteristic parameters.
- The number of concentrated masses: n=11;
- The inertial moments of concentrated masses (kgm2):
Ii = [9181.0 7999.0 7999.0 7999.0 7999.0 7999.0 7999.0 5422.0 9538.3 857.9 48240.1];
- The stiffness of shafts (MNm/rad)
Ki=[1201.92 940.82 940.82 940.82 940.82 940.82 1434.93 934.24 77.16 101.04];
- Diameters of shafts (mm):
Di =[670 670 670 670 670 670 670 670 500 570];
Designating unknown damping coefficients ai and bi,i+1, with symbles s, z, f, y, r and place them into
equation (1), we get matrix FK(22,22) for each value of kw (k- harmonic order, w- rotation speed of
propeller), having form:
Fk =
[1201920-9.18*kw2, -1201920, 0, -kw*r, kw*r, 0]
[ -1201920,2142740-7.99*kw2, -940820, 0,0,0,0,0,0,0,0, kw*r, -kw*z, kw*s, 0]
[ 0, -940820,1881640-7.99*kw2,-940820,  kw*s,-f*kw, kw*s, 0]
[ -kw*r, kw*z, -kw*s,.. ,-1201920, 2142740-7.99*kw2,-940820, 0]
[ 0, 0, ... ,-77160, 178200-0.85*kw2, -101040]
[ 0,  0,.. ,kw*y, ... 0, -101040, 101040-4.82*kw2]
To determine the determinant det(Fk) and the algebraic complements Di,j can be done with using the
Symbolic Math Toolbox package of Matlab software by the command “det(Fk)” if the matrix size is
small enough (depending on the number of non-zero elements), the result is a formula in symbolic form
on the screen, for example:
ans=
2871048086689088639517089901.315*kw14*r3*s4*z+
2218063572479354423019604109325.3*kw16*r2*s3*z +  +12883815368824962090418376
Finally, with each value of kw, after having computed all the algebraic complements Di,j of elements in
the interesting lines, replacing them into equation (8), supposing we get an equation as:
15 2 7 4
1 2 1 0... ....i p p pC f r z C z s C s C y C r C . (9)
  1,1...2i n  ,1...2i n  
 , ,i j j iD ,j iD
 1 2 1 21 2 1 2. ... . ...
m m n n
ix a a b b
416
There are p constants Ci in the left side of the equation (9), according to k and rotation speed w of the
propeller shaft. Here, to determine r, f, z, we should obtain p values of torsion torque amplitudes
according to p value of kw in order to make p equations (9) and a system of p linear equations with p
hidden: x1=r, x2=f, x3=z, , xi=fm rnz, Resolving this system we would identify the hidden r, f, s
If the computing is performed correctly and with inputs being the elastic torque values obtained by
calculations as described above (assuming the damping coefficients are constant), then the result must
match, it means, for example:
if: x1=r, x2=f, x3=z, xi=fm rnz,
then: xi=x2mx1nx3.
Being sure it correct, take values of elastic torsional torques measured practically to put into (9) to
calculate. If xi≠x2m x1nx3, for example, depending on differences we may assess the suitability of the used
assumptions about resistances.
However, in the example case, unfortunately we can not directly use the command det(Fk). It would
cause errors “Out of Memory” or the text is too long, we would get a report on screens, for example:
ans=
2871048086689088639517089901.315*kw^14*r^3*s^4*z+
  +.2218063572479354423019604109325.3*kw^16*r^2*s^3*z +
12883815368824962090418376... Output truncated. Text exceeds maximum line length of 25.000
characters for Command Window display.
To overcome this problem, we should divide the calculation procedure into some steps:
- Step 1: replace the original terms of the first matrix by new symbols (to get more compact form),
for example:
I1=-I1*kw^2+K1; r=kw*r; I2=-Ic*kw^2+K1+Kx; z=kw*z; s=kw*s; I3=-Ic*kw^2+2*Kx;
I4=-Ic*kw^2+2*Kx; I5=-Ic*kw^2+K3+Kx; I6=-I3*kw^2+K3+K4=I6; I7=-
I4*kw^2+K4+K5; I8=-I5*kw^2+K5+K6; I9= -Ic*kw^2+K6,
to obtain matrix Fk looking as:
Fk=
[ I1, -K1, 0, 0,0,0,0,0,0,0,0,-r, r,0,0,0,0,0,0,0,0,0];
[-K1, I2,-Kx, 0,0,0,0,0,0,0,0, r,-z, s,0,0,0,0,0,0,0,0];
..
[r,-r,0,0,0,0,0,0,0,0,0,I1,-K1,0,0,0,0,0,0,0,0,0];
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-K5, I8,-K6];
[0,0,0,0,0,0,0,0,0,0,y,0,0,0,0,0,0,0,0, 0,-K6, I9];
- Step 2: use the comand “det” for a smaller size, for example: D10=det(Fk(1:10,1:10)). The
obtained results looks like that:
D10 =
- I3*I4^3*K1^2 - I1*I4^3*Kx^2 + I1*I2*I3*I4^3 + I4^2*K1^2*Kx^2 -
I1*I2*I4^2*Kx^2 + 2*I3*I4*K1^2*Kx^2 + 2*I1*I4*Kx^4 - 2*I1*I2*I3*I4*Kx^2 -
K1^2*Kx^4 + I1*I2*Kx^4
- Step 3: transform Fk to Fk1 by Gause method so that Fk1(1:10,1:10) be a triangle matrix;
417
- Step 4: replace elements of Fk1(11:22,11:22) with an other symbols to get a new matrix of
12×12:
Bk=
[ b1_1, b1_2, b1_3, b1_4, b1_5, b1_6, b1_7, b1_8, b1_9, 0, 0, b1_12]
[ b2_1, b2_2, b2_3, b2_4, b2_5, b2_6, b2_7, b2_8, b2_9, 0, 0, 0]
[ b3_1, b3_2, b3_3, b3_4, b3_5, b3_6, b3_7, b3_8, b3_9, 0, 0, 0]
[ b4_1, b4_2, b4_3, b4_4, b4_5, b4_6, b4_7, b4_8, b4_9, 0, 0, 0]
[ b5_1, b5_5, b5_3, b5_4, b5_5, b5_6, b5_7, b5_8, b5_9, 0, 0, 0]
[ b6_1, b6_2, b6_3, b6_4, b6_5, b6_6, b6_7, b6_8, b6_9, 0, 0, 0]
[ b7_1, b7_2, b7_3, b7_4, b7_5, b7_6, b7_7, b7_8, b7_9, 0, 0, 0]
[ b8_1, b8_2, b8_3, b8_4, b8_5, b8_6, b8_7, b8_8, b8_9, 0, 0, 0]
[ b9_1, b9_2, b9_3, b9_4, b9_5, b9_6, b9_7, b9_8, b9_9, b9_10, 0, 0]
[ 0, 0, 0, 0, 0, 0, 0, 0, b10_9, b10_10, b10_11, 0]
[ 0,  0, 0, 0,  . b11_10, b11_11, b11_12]
[b12_1, 0, 0, 0, 0, 0, 0, 0, 0, 0, b12_11, b12_12]
- Step 5: repeat step 1), for example: det(Bk(1:6,1:6)), and then do the same as in step 2)
The purpose of this algorithm is to divide the final formula of det(Fk) into shorter terms so that they can
be displayed on the screen.
3. Conclusion
Thus, the above presented algorithms allow determining the damping coefficients if having measured
torsion elastic torque amplitudes at any shaft section. Numbers of measurement results should obtain
will depend on number of masses.
Now the paper’s author has just got the formula for det(Fk) in a form of string of characters written with
syntax errors for Matlab, which is to be corrected and the final results haven’t obtained yet because of
that it takes long time, but hopping reach it in the near future. Once, when we obtain this formula already
in an analytic form, we may use it for other shaftings having the same or less numbers of concentrated
masses.
References
[1] Đặng Hộ, Thiết kế trang trí động lực tàu thủy (tập 2), Vietnam Transport Publisher, 1986.
[2] Documentation of torsional vibration computations and measurements submitted to the GL for
approval by Polish ship design Agency, Propulsion System Torsional Vibration Analysis B 170-V.
[3] Nguyễn Mạnh Thường, Tính nghiệm dao động xoắn cho hệ trục loạt tàu B 170 –V, Vietnam
Transport Publisher, 4-2011.
[4] Nguyễn Mạnh Thường, Nghiên cứu xây dựng phương pháp xác định các hệ số sức cản bằng thực
nghiệm phục vụ tính toán dao động hệ trục tàu thủy. Scientific seminar on Shipbuilding, Water
Transport, Marine oil Industry& Equipment, Road& Rail Way Transport, Transport Publisher,
Hanoi, 4-2014, p. 34-46.