Inverse modeling approach for transformation of propeller shaft angular deformation and velocity to propeller torque load

Marine Structures 67 (2019) 102614

Contents lists available at ScienceDirect

Marine Structures journal homepage: www.elsevier.com/locate/marstruc

Inverse modeling approach for transformation of propeller shaft angular deformation and velocity to propeller torque load

T

Dražen Polića,b,∗, Sören Ehlersc, Vilmar Æsøyb a b c

Department of Ocean Operations and Civil Engineering, Norwegian University of Science and Technology, Ålesund, 6025, Norway Department of Marine Technology, Norwegian University of Science and Technology, Trondheim, 7491, Norway Institute of Ship Structural Design and Analysis, Hamburg University of Technology, Hamburg, 21073, Germany

ARTICLE INFO

ABSTRACT

Keywords: Ice-related torque load Propeller shaft response Inverse problem Transmissibility of displacement and force Finite-mode discretization Artificial noise

In full-scale trials or model-scale tests, the propeller torque load is typically measured indirectly on the propeller shaft. Due to structural flexibility and mass of the propeller and propeller shaft, the measured propeller shaft torque response includes a time delay and amplitude difference compare to the propeller load. These differences are usually small during propeller operation, but may be significant during transient states (i. e. starting of the engine, rough sea, or operation in ice-covered waters). In this paper, an inverse model is used to calculate propeller load based on the propeller shaft response. The results indicate that the inverse model can be less complex than the complete propulsion machinery model and consider only the propeller and propeller shaft. Several modeling approaches are tested based on four proposed measuring setups and three forms of system discretization. Out of five inverse modeling approaches, two approaches showed the ability to transform measured propeller shaft angular deformation and velocity responses, colored with noise, into the corresponding rule-based propeller torque load with acceptable accuracy. Furthermore, the robustness of the inverse approaches to the random ice-propeller interaction, sampling frequency, measuring location along the propeller shaft and starting time (before and during ice-propeller interaction), is discussed.

1. Introduction The most common source of ship propulsion is a propeller, which is traditionally optimized for operation in calm sea conditions with a steady load [1]. However, during ship navigation in waves or ice-covered waters, large transient propeller loads can occur. Ice impact on the propeller could be the most extreme case of such loads. At present, our understanding of ice-related transient loading mechanisms, based on measurements and simulations, is not necessarily sufficient for informed and safe design of the propeller and other elements in the mechanical system (e. g. propeller shaft, diesel engine or electrical motor). The requirements, advantages and challenges for each method used for measuring and simulating propeller loads in the presence of ice are listed in Table 1. The applicability and fidelity of each method, given in Table 1, are estimated and illustrated in Fig. 1. The applicability of propeller load measurements depends on the information collected during the full-scale trial or model-scale test: ice properties,

∗ Corresponding author. Department of Ocean Operations and Civil Engineering, Norwegian University of Science and Technology, Ålesund, 6025, Norway., E-mail address: [email protected] (D. Polić). URL: https://www.ntnu.edu/sast (D. Polić).

https://doi.org/10.1016/j.marstruc.2019.04.002 Received 28 September 2018; Accepted 8 April 2019 0951-8339/ © 2019 Elsevier Ltd. All rights reserved.

Marine Structures 67 (2019) 102614

D. Polić, et al.

Table 1 Methods for prediction of the propeller load in the presence of ice. Options

Requirements

Advantage

Challenges/Disadvantages

Full-scale propeller load measurements (e.g. Refs. [2–10])

- A specially-designed propeller with holes for the strain gauges and instrumentation cables. - A ship instrumented with measuring devices for measuring ice-propeller interaction (e. g. ice-properties, environmental and operational conditions).

- First-hand information.

- Direct exposure of the strain gauges, attached to the propeller blade surface, to ice impacts. - Maintenance access in the field is limited. - High installation costs due to additional grooves in the propeller and propeller shaft for the instrumentation cables. - Generally ice loads acting on the propeller cannot be evaluated in terms of ice properties, operational and environmental conditions [11]. - Rare measurement setup.

Model-scale propeller load measurements (e.g. Refs. [12–15])

- An ice tank. - A self-propelled model-scale ship. - Blade dynamometer.

- Ice properties, and operational and environmental conditions are known. - Six component blade dynamometer can be used for measuring propeller blade load [13,14,16].

- The internal mechanics of model-scale ice differ from sea ice [17,18] - The prime mover of the propeller in model-scale tests is an electrical motor, giving no information about a full-scale coupled relationship between the diesel engine and directly driven propeller used by most merchant ships. - Hydrodynamic scaling effects.

Inverse propulsion machinery (PM) models [19,20]

- Full-scale or model-scale propeller shaft measurements.

- Transforms reliable, cheap and common propeller shaft shear strain (angular deformation) measurements to the propeller load.

- Full-scale trial or model-scale test is required. - The required regularization method needs additional information about the PM system, which has to be assumed. - Computationally expensive. - The sensitivity of the inverse PM model to uncertainty (e. g. noise) in the measured signal is unknown. - Inverse PM models [19,20] are not validated against model-scale or fullscale simultaneously measured propeller load and propeller shaft response. - Only the torsional component of the propeller load is transformed.

Ice-propeller prediction models (e.g. Refs. [11,21–26])

- Easily applicable empirical models of ice failure, usually derived from specialized tests.

- Universal application (wide range of propeller designs and ice properties). - Inexpensive (no need for a fullscale trial or model-scale test).

- Validity of the method depends on empirical factors and weight functions derived from simple and fundamental tests that oversimplify ice-propeller interaction. - Unsteady ice and propeller motion is simplified: ○ an ice shape is simplified(e. g. sphere, infinitive plate). ○ propeller is considered as rigid body with constant angular and translational velocity. ○ propulsion machinery response is ignored. ○ interaction between ice block and water is simplified or neglected. - Validated against limited propeller load data collected through full-scale trials and model-scale tests. - Often computationally expensive.

operational (e. g. RPM, ship speed and bearing) and environmental (e. g. ice block size and shape, trajectory and number of ice blocks in front of the propeller) conditions. The fidelity of the measured propeller load depends on the design of the measurement setup, data loggers, sensors, noise filters, and in the case of model-scale tests, validity of the model-scale ice. Ice-propeller prediction models, which nowadays set the design load in the rules, are universally applicable because ice properties, as well as operational and environmental conditions, are inputs to the model. The fidelity of ice-propeller prediction models depends on the propeller load measurements that are used for validation, and the validity of empirical factors and weight functions derived from specialized tests. 2

Marine Structures 67 (2019) 102614

D. Polić, et al.

Fig. 1. Assessment of methods used for prediction of the propeller load in the presence of ice. Case- and trial-specific application represents the level of information collected about ice properties, operational and environmental conditions. For trial-specific application, the aforementioned is always incomplete.

The propeller load calculated using inverse propulsion machinery (PM) models depends on full-scale trials or model-scale test propeller shaft measurements, so the applicability is the same as for propeller load measurements. Compared to the propeller load measurements, the fidelity of propeller load calculated using inverse PM models is much lower. The reason for low fidelity is not the reliability of the propeller shaft measurements, but the fact that the inverse PM models developed by Refs. [19,20] are not validated against model-scale or full-scale simultaneously measured propeller load and propeller shaft response (angular deformation and velocity). Furthermore, the sensitivity and stability of the inverse PM models with respect to error (noise) in the measured signal is not clearly presented in Refs. [19,20]. 1.1. Propeller shaft measurements history and state of the art inverse PM models The history of using propeller shaft measurements (see e. g. Refs. [7,20,27–32]) to estimate propeller load in the presence of ice dates back to the 1960's [33]. The ice regulations given by Finnish-Swedish Ice Class Rules (FSICR) and Det Norske Veritas (DNV) between 1971 and 2008 for propeller design were based on corrected/transformed propeller shaft measurements. Other than [19,20], there is, however, little available information regarding how to transform these measurements and which assumptions or simplifications to employ. Browne et al. [19] invert the measured propeller shaft torque Qs to obtain the propeller torque (Qp ) load using a deconvolution method, a technique widely used in signal and image processing, without using any simplifications or assumptions. The shape of the calculated Qp load (amplitude, duration, and timing of individual blade loads) is obtained by determining the series of impulses, which are summed up in the measured Qs . Ikonen et al. [20] present an alternative inverse propulsion machinery (PM) model. The dynamic behaviour of the system in Ref. [20] is modeled according to the idea presented by Ref. [34], where an incremental form of the governing equation of torsional vibration of the shaft is combined with the Newmark-β time integration method. The system of governing equations is inverted with help from a regularization method. Regularization methods can also be combined with the previously described de-convolution method, because the result of the inversion of impact response may not always be well-posed (unique, stable and existing), according to Refs. [35–37]. The disadvantage of the regularization methods and hence inverse models is that the required additional information is often derived from (arbitrary) assumed output and, consequently, the transformation procedure is only valid if the result is identical or similar to the assumed output. An additional challenge in both inverse procedures is the inversion of large matrices since the Qp at each time step is calculated based on roughly 500 measured data points. 1.2. Objective and hypothesis In this paper, a new inverse model is developed. The governing differential equation of torsional vibration is simplified to a solvable system of equations that do not require a regularization method, if the measured time histories of propeller shaft angular deformation ( s ) and velocity ( s ) are continuously differentiable in time. The number of data points that have to be considered in the calculation of the propeller torque load at a certain time instant is reduced from 500 to 2 points, and the inverse model is reduced by considering only the propeller and propeller shaft (two sub-models) instead of the whole PM system illustrated in Fig. 2. The proposed inverse model is numerically validated by using a separate model called the PM model. The PM model is adopted from Ref. [38] and it considers all torsional mechanical elements in the directly driven PM system. The PM model is used for calculation of the 3

Marine Structures 67 (2019) 102614

D. Polić, et al.

Fig. 2. Propulsion machinery (PM) with four different setups for propeller shaft response measurements.

propeller shaft angular deformation s and velocity s (which are used as inputs in the inverse model) based on the rule-based and random propeller torque. To avoid confusion between propeller torque used by the PM model as input and propeller torque calculated by the inverse model, the calculated propeller torque is called inverse propeller torque and designated Qp* . The inverse torque Qp* is compared to the torque Qp , the input to the PM model, and the residual and the number of oscillation cycles within both torques are used to quantify the accuracy. Furthermore, the sensitivity of the residual to: -

Different torque Qp inputs Artificial noise (white, pink, red) that is added to the propeller shaft angular deformation The four different input sets (propeller shaft measurement setups) illustrated in Fig. 2 The four different sampling frequencies Two different starting of transform is presented.

s

and velocity

s

1.3. Outline The inverse procedure and inverse model are described in Section 2, followed by the validation procedure and PM model in Section 3, and then a sensitivity study in Section 4. Finally, the most suitable measurement setup, based on the sensitivity analysis, is discussed in Section 5, together with a future development and possible application. 2. Inverse methodology The basic concept of the inverse methodology is to calculate the causal impact, the Qp* load, that produced the observed propeller shaft angular deformation s and velocity s signal data. Fig. 3 illustrates such an inverse process, which comprises three steps: - Processing (filtering) of the measured signals that include noise - Correction of the filtered angular deformation (torque) and velocity for the propeller shaft dynamics (inverse propeller shaft submodel) - Correction of the calculated propeller shaft torque at the end connected to the propeller for the propeller dynamics (inverse propeller sub-model). 2.1. Measurement setups The measured data must be acquired from accessible location xk along the propeller shaft. Since model-scale or full-scale data were unavailable, simulation based dataset is calculated from the PM model presented in Section 3, colored with the artificial noise, and used to test and verify the inverse methodology. Four types of measured data are assumed as possible inputs. The first option is a time history record of angular deformation s (xk , t ) , the second option is the deformation s (xk , t ) and angular velocity s (xk , t ), and the last two options are set by collecting an additional angular deformation s (xk , t ) signal compared to the first and second option, respectively. In each option, the number of measuring locations xk along the propeller shaft must be proportional to the number of signals. The four assumed measured data sets are provided by measurement setups 1–4, illustrated in Fig. 2. Typical measurements include a single propeller shaft torque derived from angular deformation (see e. g. Refs. [7,20,27–32]) and in some cases a lowfrequency angular velocity (see e. g. Refs. [20,29,31,32]). 4

Marine Structures 67 (2019) 102614

D. Polić, et al.

Fig. 3. Inverse methodology.

2.2. Filtering the noise Processing of the measured data starts with acknowledging uncertainty in the measurement (spurious readings, measurement error, background noise, etc.). The interfering noise can be partially or completely suppressed by a filter, designed based on practical and/or numerical knowledge about the measured response and possible noise. Typically, a low-pass filter (LPF) would be applied to keep the low-frequency signal of interest and eliminate noise above a certain cutoff frequency. Here, due to a lack of specific knowledge about the noise, the LPF is replaced with a more general locally weighted regression [39], that smooths the signal over a wide range of frequencies. The effect of regression on the measured angular deformation (torque) is illustrated in Fig. 4. After smoothing (filtering), the measured data is used as input to the inverse model. 2.3. Inverse model Excluding the possibility of an off-design plastic deformation and assuming a uniform propeller shaft cross-section, a measured angular deformation s (xk , t ) can be expressed as a torque

Qs (xk , t ) =

2Ip G Dshaft

s (x k ,

t ).

(1)

The constants Ip , Dshaft and G are the propeller shaft cross-section second polar moment of area, outer diameter, and shear modulus, respectively. According to Ref. [40], the torque Qs (xk , t ) differs from Qp (t ) due to the propeller and propeller shaft dynamics

Fig. 4. Influence of white noise and regression on the simulated torque Qs (xk , t ) during steady state condition before ice-propeller interaction. 5

Marine Structures 67 (2019) 102614

D. Polić, et al.

Fig. 5. Difference between propeller torque Qp (t ) load and propeller shaft torque Qs (xk , t ) response at locations xk .

as illustrated in Fig. 5. Hence, the purpose of the inverse model is to reduce the difference by accounting for the dynamic influence of the propeller shaft (inverse propeller shaft sub-model) and propeller (inverse propeller sub-model) to a certain extent.Table 2 shows five possible approaches (ai) with respect to four measurement setups (illustrated in Fig. 2). Three approaches (a2, a3 and a4) consider the complete contribution from propeller and propeller shaft dynamics, one approach (a1) consider a partial contribution (only dynamic influence of the propeller), and one approach (a0), the base case, does not considers any dynamic contribution by equalising signal Qs (xk , t ) with Qp* (t ) . As shown in Table 2, four approaches are redundant because the full potential of the measurement setup is not utilised and three approaches are infeasible: the time history of angular velocity s (xk , t ), required by the inverse propeller submodel, is not measured in setup 1 and 3; and number of measured signals is less than the number required by the inverse propeller shaft sub-model.In the case of approaches a1, the inverse propeller shaft sub-model is neglected and the measured data is directly used as an input to the inverse propeller sub-model. A partial contribution of the propeller shaft dynamics, opposite to the approach a1, is not considered in the paper because the torque variation Qs (xk , t ) over the shaft length is relatively small compare to the propeller contribution. Furthermore, approaches a1 and a2 are improved and corrected versions of the rigid and flexible inverse model presented by Ref. [41] (Note: inverse model in Ref. [41] is referred as a reverse model). 2.3.1. Inverse propeller shaft sub-model The inverse propeller shaft sub-model is the first distinct part of the inverse model and calculates torque Qs (t ) and angular velocity s (t ) at the boundaries x i = [0, Lshaft ] based on the smoothed measured data collected at x k . Hence, the equation of torsional vibration

Js ¨s (x , t ) + cs s (x , t ) + ks s (x , t ) =

(2)

Qs,2 (t ) + Qs,1 (t )

should be re-organized according to the measurement setup (data inputs s and s ). Indices s,1 and s, 2 are donated to the torque load Qs (t ) at the boundary x i = Lshaft (the end connected to the flexible coupling) and x i = 0 (the end connected to the propeller), respectively. The propeller shaft and its connections are illustrated in Fig. 2. The constants Js , cs and ks in Eq. (2) are propeller shaft inertia, damping, and stiffness, respectively.Re-organization of Eq. (2) starts by deriving the governing differential equation of torsional vibration as a system of differential equations by considering a finite-mode method. In contrast, the inverse PM models presented by Refs. [19,20] consider the finite-lump and finite element methods. The finite-mode method uses a separation of variables (Fourier method) to discretize the shaft into a series of imaginary modal shafts that are connected in a parallel “circle”. This type of harmonic decomposition, an alternative approach to the spatial distribution that is incorporated in the finite element and finite-lump methods, yields one of the most accurate low order models used for the simulation of wave-like response [42].For forced or unforced beam response, the propeller shaft angular position ( s ) can be expressed as a sum of kinematically valid orthogonal mode shapes (Yn ) multiplied by unknown time-dependent generalized coordinates ( n ) [43–45]: s (x ,

t) =

Yn (x )

n (t )

(3)

n=0

Substituting Eq. (3) into Eq. (2), multiplying each term in Eq. (2) by the mth mode shape (Ym ), and integrating over the shaft Table 2 Approach (a) used in the inverse model and dynamic consideration that are taken into account (✓). Infeasible approaches are mark with minus symbol (−), while redundant approaches are mark with plus symbol (+). Shaft dynamics Propeller dynamics Measurement Measurement Measurement Measurement

setup setup setup setup

1 2 3 4

✗ ✗

✗ ✓

✓ ✓

a0 + + +

– a1 – +

– a2 a3 a4

6

Marine Structures 67 (2019) 102614

D. Polić, et al.

length (Lshaft ), we obtain

(

)

Lshaft

Jn (x )dx ¨n (t ) +

0

(

Lshaft 0

cn (x )dx

)

n (t )

+

(

Lshaft 0

)

kn (x )dx

n (t )

=

Qs,1 (t )

Lshaft 0

(x ) Ym (x )dx + Qs,2 (t )

Lshaft 0

(x ) Ym (x )dx,

(4) where the only nonzero contributions, as explained below, are for n = m . The right hand side of Eq. (4) arises from the expression for an external point torque such as Qs (t ) . Torque Qs acting at location x i can be expressed using the delta function (x ) [46]:

0, x xi 1, x = [0, Lshaft ],

(x ) =

Lshaft

and hence the integral 0 orthogonality of the modes Lshaft

(x ) Ym (x ) dx is equal to the mode shape value Ym (x ) at location x = [0, Lshaft ]. Accounting for the

Yn (x ) Ym (x ) dx = 0 for n

0

(5)

m,

(6)

the modal inertia ( Jn ) is equal to Lshaft

Jn =

0

Ip

Yn (x ) Ym (x )dx = n=0

Lshaft 0

Ip Yn2 (x )dx,

(7)

and the modal stiffness (kn ) is equal to 2 n

kn =

L shaft 0

Jn (x )dx.

(8)

The constant ρ is the density, while n is the modal natural frequency. The second polar moment of area Ip is not a function of x because the cross-section of the shaft is assumed to be constant over the shaft length. The modal damping (cn ) is omitted in the rest of the paper because damping has negligible dynamic contributions far from resonance, as expected for a steel shaft under the given loads. Thus, the equation of motion considering modal decomposition can be written as

Jn ¨n (t ) + kn

n (t )

=

(9)

Qs,2 (t ) Yn (0) + Qs,1 (t ) Yn (Lshaft ).

If the right side of Eq. (9) can be decomposed into the sum of simple functions, i. e. i j Qi, j cos( i, j t ) Yn (xi ) where Qi, j represents the amplitude of the external torque and i, j denotes its frequency, the importance of modal inertia and stiffness in a certain mode can be checked using the equation of transmissibility of vibration 2 i, j 2 n

2 i, j

+

2 n

2 n

2 i, j

= 1.

(10)

Indices i and j represent the number of external torques at the boundary (in this case Qs,1 and Qs,2 ) and the number of harmonic components, respectively. Eq. (10) is obtained by dividing Eq. (9) with the right hand side. The first term in Eq. (10) is called transmissibility of force (TRF) and describes a dimensionless measure of how much of the external torque of amplitude Qi, j is transmitted from the point of the application (boundary) to the inertia Jn . The second term in Eq. (10), called transmissibility of displacement (TRD), describes how much of the external torque of amplitude Qi, j is transmitted to the displacement n . Fig. 6 illustrates TRF and TRD for an undamped dimensionless shaft as a function of the frequency ratio i, j / n . Near i, j / n = 1, or resonance, the maximum amount of the external torque is transmitted to the inertia and displacement, respectively. Assuming small frequency ratios i, j / n 0.2 as expected for the PM system, transmission of external torque to displacement will be at least 25 times larger than the transmission to inertia. Hence, for flexible modes, the equation of motion (Eq. (9)) can be simplified to.

Fig. 6. A comparison between TRF and TRD for undamped dimensionless shaft on a semi-logarithmic plot. 7

Marine Structures 67 (2019) 102614

D. Polić, et al.

Qs,2 (t ) Yn (0) + Qs,1 (t ) Yn (Lshaft )

n (t )

for n > 0,

kn

(11)

or in a more general form n (t )

Qs, i (t ) i

Yn (x i ) for n > 0. kn

(12)

The location x i represents points along the shaft where external torque Qs is applied. Furthermore, by knowing that the rigid mode does not contribute to the measured angular deformation location xk multiplied by the radius), the deformation s can be expressed as s (x k ,

t) =

Dshaft

s (x k ,

2

t)

x

=

Dshaft 2

Y n (x k )

t)

Qs,2 (t )

Dshaft 2

n>0

(a twist at a certain

n (t ).

(13)

n>0

Substitution of the simplified solution of displacement s (x k ,

s

n

(Eq. (11)) into the equation of angular deformation (Eq. (13)) yields

Dshaft Y n (xk ) Yn (0) + Qs,1 (t ) kn 2

Y n (xk ) Yn (Lshaft ) n> 0

kn

,

(14)

or in a more general form s (x k ,

Dshaft

t)

2

Qs, i (t ) i

n> 0

Y n (xk ) Yn (x i ) . kn

(15)

Y n (x ) is the first spatial derivative of mode shape Yn (x ) .The negligible TRF contribution also implies that Jn ¨n (t ) 0 in the flexible 0 , ¨n (t ) 0 . Hence, the second derivative of the right hand side of Eq. (3) can be simplified to one modes (n > 0 ), and because Jn term: ¨s (x , t )

(16)

Y0 (x ) ¨0 (t ).

Substituting ¨0 (t ) from Eq. (9) for the rigid mode into Eq. (16), the first derivative of the measured angular velocity be expressed as s (x k ,

t)

Y0 (xk )

t

s (x k ,

t ) , can

Qs,2 (t ) Y0 (0) + Qs,1 (t ) Y0 (Lshaft ) (17)

J0

or in more general form s (x k ,

t)

Y0 (xk ) J0

t

Qs, i (t ) Y0 (x i ).

(18)

i

From a system of independent linear equations, represented by Eq. (15) and Eq. (18), i unknown torques Qs, i can be found if the number of equations is equal to k, the number of measured variables ( s and s ). Therefore, for finding torques Qs,1 and Qs,2 that act at the propeller shaft ends, two time history signals must be measured (e.g. Setup 2 and 3). The third time history signal, propeller shaft angular velocity s (x3 , t ) measured in the setup 4, is not used in the system of equations represented by Eq. (15) and Eq. (18), but in a control of propeller shaft angular acceleration ¨s,2 (t ) output from the inverse propeller shaft sub-model. The index s, 2 of angular acceleration ¨s, i (t ) indicates the same boundary as in the case of torque Qs,2 , and control of the ¨s, i (t ) is explained later. After solving Eq. (15) and Eq. (18), solution Qs, i is used further for finding the propeller shaft angular acceleration s, i / t at location x i . The angular acceleration s, i / t is calculated from the second time derivative of Eq. (3), where the contribution of flexible modes (n > 0 ) to the angular acceleration is found from the second time derivative of Eq. (12) and the contribution of the rigid mode to the angular acceleration is found from Eq. (9) for n = 0 . Thus, the equation for angular acceleration s, i / t can be express as

t

s , i (t )

Y0 (x i ) ¨0 (t ) +

Yn (xi ) n>0

1 2 kn t 2

Qs, i (t ) Yn (x i ),

(19)

i

where

¨0 (t ) =

1 J0

Qs, i (t ) Y0 (x i ).

(20)

i

In setup 4, where one of the measured k variables is s (xk , t ) and not used for calculation of Qs, i torque, the contribution of flexible modes (n > 0 ) is subtracted from the first time derivative of s (xk , t ) . Hence, a new set of ¨0 (t ) data can be derived from

¨0 (t ) =

1 Y0 (xk )

t

s (x k ,

t)

Yn (xk ) n>0

1 2 kn t 2

Qs, i (t ) Yn (xi ) , i

(21)

and used to control the previous set of data ¨0 (t ) derived from Eq. (20). The contribution of flexible modes is found as earlier, from the second time derivative of Eq. (12).Compare to the previous equations (e.g. Eqs. (15), (18) and (19)) that may consider interpolated values of the measured velocity s (xk , t ) (e.g in the case, when all measured signals are not sampled with the same sampling 8

Marine Structures 67 (2019) 102614

D. Polić, et al.

Fig. 7. Linear correction and control of ¨0 (t ) using derived from Eq. (21).

s (xk ,

t ) in a4. Modal angular acceleration ¨0 (Qs, i) is derived from Eq. (20), while ¨0 (Qs, i,

s)

is

frequency), Eq. (21) considers only time instances that are measured. Hence, modal angular acceleration ¨0 (t ) derived from Eq. (9) for n = 0 in Eq. (20) is linearly interpolated to go through the points derived by Eq. (21). Use of such interpolation can be justified by the fact that sensitivity of the solution ¨0 (t ) to potential noise is smaller in Eq. (21) than in Eq. (20) because Qs, i (t ) is divided with the stiffness kn that is greater than J0 . Furthermore, considering that propeller shaft torque Qs (xk , t ) is typically measured with a higher sampling frequency than angular velocity s (xk , t ) , and that Qs, i (t ) can be found from Qs (xk , t ) solely, set of ¨0 (t ) from Eq. (21) alone would be too small. The principle of such control is illustrated in Fig. 7, and the benefits of the control, in the presence of noise in the torque Qs (xk , t ) and angular velocity s (xk , t ) , are presented in Section 4. In a2, Eq (21) can be used but the benefit of control of ¨0 (t ) is neglectable duo to the fact that Eq. (21) has been simplified as Eq. (18) and used for finding torque Qs, i . 2.3.2. Inverse propeller sub-model The second distinct part of the inverse model is the propeller sub-model, which calculates the torque Qp* (t ) based on the output from the inverse propeller shaft sub-model: the torque Qs,2 (t ) and angular acceleration ¨s,2 (t ) . Unlike the shaft sub-model, the inverse propeller sub-model considers only rigid body discretization because the propeller is relatively stiff in torsion compared to the shaft. The dynamic contribution from inertia is calculated by multiplying the propeller moment of inertia (Jp ) with the first derivative of s,2 (t ) , and the obtained product is subtracted from the torque Qs,2 (t ) , giving:

Qp* (t ) = Qs,2 (t )

Jp

t

s,2 (t ).

(22)

3. PM and validation methodology The inverse methodology illustrated in Fig. 3 is numerically validated by using an additional model, called PM model, which calculates the propeller shaft response (input to the inverse model) based on a known propeller torque load Qp . The calculated propeller shaft response ( s and s ) is disturbed by artificial noise and applied to the inverse model (a0-a4). The known torque Qp is compared with the torque Qp* output from the inverse model, and the residual together with the number of half cycles obtained using rain-flow counting algorithm is used as a criterion of the accuracy of the inverse methodology. This validation process is illustrated in Fig. 8. 3.1. PM model Unlike to the inverse model, the PM model considers all components of the PM system: the fixed pitch (FP) propeller, propeller shaft, flexible coupling and two-stroke slow-speed seven-cylinder diesel engine. In the PM model adopted from Ref. [38] and illustrated in Fig. 9, the PM components and their connections (illustrated in Fig. 2) are formalized as a network of idealized processes (e. g. energy supply, storage, dissipation, conversion or routing) using the bond graph (BG) method [46–48]. Other than the propeller shaft, which is discretized using the finite-mode method for higher accuracy, all PM components are spatially distributed using a traditional finite-lump method. Inertias, stiffnesses and damping are represented with I-, C- and R-elements, respectively. For more information about representation of the finite-lump and finite-mode methods using BG-elements (I-, C- and R-) see Refs. [46–50]. Application of the BG method and block digram method (bi-directional exchange of energy used by BG method is replaced with unidirectional flow) for transient simulation of the PM system response can be also found in Refs. [51–55]. In the PM model, the FP propeller is considered as an equivalent rigid disc which is connected to the diesel engine through the propeller shaft and a linear flexible coupling with a static stiffness and damping. Assuming that the propeller and flexible coupling are not connected at the very edge of the propeller shaft, free-free boundary conditions with angular deformation equal to zero at the both ends are considered. The mode shape function and modal natural frequency for free-free boundary conditions are: 9

Marine Structures 67 (2019) 102614

D. Polić, et al.

Fig. 8. Validation of the inverse methodology. Note: propulsion machinery (PM) model is presented in Ref. [38].

Yn (x ) = cos

n x Lshaft

(23)

and n

=

n Lshaft

G

.

(24)

The cross-section of the propeller shaft is assumed constant and set by outer diameter Dshaft and inner diameter dshaft . The propeller and flexible coupling are connected to the propeller shaft at location x i = [0, Lshaft ], where Lshaft is a propeller shaft length. The flexible coupling inertia is added to the engine's flywheel, together with the rear journal inertia from crankshaft. The cylinder units are connected in-line with an equivalent crank throw spring and damper estimated as in Ref. [56], the average of Carter's and Wilson's formulas. Complete rigid body motion of the crank mechanism (piston, piston rod, crosshead and connection rod) is included in the cylinder unit through the IC-field, a special BG-element [46,47]. The cylinder pressure ( pcyl ) that acts on the top of the piston is

10

Marine Structures 67 (2019) 102614

D. Polić, et al.

Fig. 9. PM model adopted from Ref. [38].

calculated at each step by scaling the cylinder pressure measured at the engine's maximum continuous rating (MCR) with a Proportional-Integral-Derivative (PID) controller.1 After the scaled pressure curve is determined, the corresponding cylinder pressure value for each cylinder unit is estimated from the crankshaft position. The advantage of the proposed modeling technique, in the case of simulation of transient PM response during the ice-propeller interaction, over a conventional lump shaft model and an equivalent disc type crank moment of inertia2 used by Refs. [29,53,57–59] is given in Ref. [38]. It is important to note that the PID controller is not validated with full-scale measurements. A complete list of all model parameters, the measured cylinder pressure curve and detailed explanations of each sub-model of the PM model (marked with a dashed line in Fig. 9), as well as the sensitivity of the system response to the number of propeller shaft modes, is given in Ref. [38]. In the present work, only the increase in the number of propeller shaft modes due to the difference between system response (global) and internal shaft response (local) is discussed in Section 3.1.1, and propeller torque variations are given in Section 3.1.2. 3.1.1. Propeller shaft - number of modes Compare to the work done by Ref. [38], where three flexible modes were enough for calculating the principal amount of energy that is exchanged at the propeller shaft ends during a transient state, the number of the flexile modes in this paper is determined by 1 2

An alternative approach to scaling the pcyl curve would be to model a diesel engine turbocharger, please see Ref. [55] as an example.

The equivalent disc type crank moment of inertia neglects the variation of inertia due to the axial or rotational oscillations of reciprocating masses, and hence the total inertia of the PM system is constant. In a engine, depending on the design, the reciprocating mass could be crankshaft throws, connection rods, crossheads, piston rods and pistons. 11

Marine Structures 67 (2019) 102614

D. Polić, et al.

Fig. 10. Convergency of Qs (xk ) torque to the value of Qsl, i torque in a steady state.

considering internal torque response Qs (xk , t ) within the propeller shaft in a steady state. During a steady state, ship operation in calm open water condition, the torque Qs (xk , t ) (Eq. (1)) and angular velocity s (xk , t ) (Eq. (3)) are constant. This makes Eq (15) from approximate into an exact solution, and the influence of the flexible modes can be expressed as

Ip G Qs ( x k ) = l Qsl, i Qs,1

Qs, i i

n>0

Y n (xk ) Yn (x i ) , kn

(25)

where

Qsl, i =

Qs , i .

(26)

i {0 xi xk }

Considering two opposing torques Qs, i that act at two ends of the propeller shaft, and substituting kn (Eq. (8)), Yn (Eq. (23)) and Y into Eq (25) yields

Qs ( x k ) = Qsl, i

2

sgn(Qs, i ) n 1sin i

n>0

n xk n xi cos . Lshaft Lshaft

n

(27)

Fig. 10a illustrates ratio Qs (xk )/ Qsl, i as a function of the flexible modes, and a convergency to the unit value for two different locations xk along the shaft. In this paper, the ratio Qs (xk )/ Qsl, i [0.975, 1.025] is considered as an acceptable inaccuracy and number of flexible modes is increased from 2 + to 53+, see Fig. 10b. The last mode (labeled with +) in the PM model is not complete due to the additional modal compliance, which omits contribution from inertia (i. e. Fig. 9 adopted from Ref. [38] illustrates the propeller shaft model with 2 + modes). For more information about modal compliance see Refs. [42,46].

12

Marine Structures 67 (2019) 102614

D. Polić, et al.

Fig. 11. Illustration of the rule-based and a random excitation sequence.

Fig. 12. Rule-based and random correlation between ice-propeller interaction angle ( i ) and impact factor (Cq) for 4-bladed propeller. The rules have the fourth case, a variation of a case with Cq = 0.5 and i = 45o witch considers a two ice pieces that simultaneously interact with the propeller. In such case, the number of impact would be double compare to illustrated DNV case 1.

3.1.2. Propeller torque load A DNV GL Ice(1A)3 ice class FP propeller with diameter (D = 6 m ), hub external diameter (d = 1.8 m ), pitch at 0.7R (P0.7 = 4.2 m ) and maximum blade thickness at 0.7R (t 0.7 = 0.102 m ) is considered. Based on the propeller geometry and mass (added mass from water and ice is neglected), the propeller inertia Jp of 46,000 kgm2 is calculated. The hydrodynamic torque (Qhyd ) is estimated as the product of a propeller loading coefficient (assumed as 8800 kgm2) and the square of the propeller angular velocity, reduced for ice operation, according to Ref. [60], to 89 rev/min (85% of MCR). The ice-related torque (Qice ) input variations (rule-based and random load) are added on the top of the Qhyd torque, input from open water condition. The rule-based Qice torque is defined according to DNV GL excitation case 1 for an open FP propeller [60]. Twelve identical halfsinus impacts and maximum ice-related torque (Qmax) of 933 kNm are calculated for a 4-bladed Ice(1A) ice class propeller. The rulebased Qmax torque is estimated based on the propeller interaction with a rectangular ice block with the dimensions Hice x 2Hice x 3Hice [60]. The maximum design ice block thickness (Hice) considered by the Ice(1A) ice class is 1.5 m. The shape of each impact, in the DNV GL case 1, is set with the ice-propeller interaction angle ( ice ) of 90 and impact factor (Cq) of 0.75, see Fig. 11. The angle ice corresponds to the propeller blade advancement through the ice block, represented in Fig. 11 as a distance between two roots. The total ice torque Qice sequence is obtained by summing the impacts from each blade while considering the angle between the blades. The sequence is also multiplied with a linear ramp function for 270 of propeller rotation at the beginning and at the end of the sequence. Other excitation cases, different ratios between αice and Cq, given by the classification societies [60–62] are illustrated in Fig. 12, but not considered here. Compared to the rule-based excitation cases that currently define the design load of the PM system [60], the random Qice torque given in this paper, as the second input to the PM model, defines the extension of the rule-based torque to account for a random number and wide range of ice impacts. A fixed number of impacts and a fixed ratio between αice and Cq are replaced with a random number and a random ratio generated using a discrete uniform distribution. The number of random ice impacts ranges from one to twelve (the maximum number according to the rules), while interval limits for ice = [5 , 135 ] and Cq = [0,1] are assumed and given in Fig. 12. The randomly generated impacts are described as simple symmetrical triangular functions appended to the sequence by

3 Det Norsk Veritas and Germanischer Lloyd (DNV GL) ice class Ice(1A) is equal to International Association of Classification Societies (IACS) ice class PC 6 and FSICR IA.

13

Marine Structures 67 (2019) 102614

D. Polić, et al.

Fig. 13. Illustration of random ice-propeller interaction. In total 400 random ice-propeller sequences. Table 3 List of assumptions used in each inverse model approach. Assumptions

a0

a1

a2, a3 and a4

cn 0 ¨n > 0 (t )

✓ ✓

✓ ✓

✓ ✓

✓

✓

✗

✗

✗

0

✓

J0 ¨0 (t )

0

kn

0

n (t )

Jp ¨s,2 (t )

0

✓

✓

✗

considering the angle between the blades, as with the rule-based excitation sequence. Fig. 11b illustrates a random ice-propeller sequence (grey line) composed of twelve independent random ice impacts (black line). In total, 400 different random ice-propeller sequences are generated and connected together “into a chain” by considering from one to five full propeller revolutions without ice-propeller interaction (Qp = Qhyd ) between each sequence, see Fig. 13. As a consequence, the random ice-propeller interaction is defined for approximately 1850 propeller revolutions or 25 min with 2623 random ice impacts. Ice-propeller interaction angle αice and ratio Cq for each random impact are illustrated in Fig. 12. Finally, the random Qice is obtained by multiplying the random ice-propeller interaction with the rule-based Qmax torque. For comparison, real-world ice-propeller interaction can be summarized as a sequence of irregular impacts, that are occasionally high. Veitch [63] observed approximately 400 ice pieces that interacted with the propeller class FSICR IA Super,4 out of which the majority, classified by the longest visible length, were smaller than 1.5 m. Only a few pieces were 3 m or longer, and between 10 and 30 pieces interacted with the propeller per minute. These observations were based on full-scale video measurements collected during three periods that in total lasted 25 min. Use of the random torque, in addition to the rule-based torque, can be justified by the fact that validation of the inverse model based on the short periodic half-sinus torque sequence may result with unrealistically good agreement. Together, the torque Qice and Qhyd give the propeller torque Qp . In the literature, Qp is sometimes called the ice-propeller torque. 3.2. Signal noise When an inverse model is numerically validated, the concept of inverse crime [64,65] should be kept in mind. The inverse crime occurs when the same model is employed to generate, as well as to invert, simulated measured data. To avoid the inverse crime, it is crucial that the PM model which is used to obtain simulated measured data has no connection to the inverse model, or the simulated measured data have to be modified to have errors in a similar fashion as the real measured data. The inverse crime is not possible in situations when full-scale or model-scale measured data are used. This is only a problem of computational studies. The presented inverse and PM models differ: simplifications of the equation of torsional vibration for each inverse model approach are marked (✓) and listed in Table 3, while none of these assumptions are made in the PM model. Nevertheless, artificial noise is added to the simulated propeller shaft angular deformation s (xk , t ) :

ˆs (xk , t ) =

s (x k ,

(28)

t ) + q (t ) e ¯s

to obtain the artificial measured deformation ( ˆs ). The mean angular deformation value ( ¯s ) is calculated based on the response, at location xk , from the open water hydrodynamic torque Qhyd produced at 85% of MCR. Furthermore, ¯s dependency on the location xk , after accounting for the inaccuracy of the PM model with 53 + flexible modes, is too small to be worth consideration. The constant e is an assumed fraction of the mean angular deformation, and q (t ) is non-dimensional noise. The non-dimensional noise is generated from uniformly distributed random numbers within the range [-1, 1] (zero mean value and standard deviation equal to 1/3 ), and transformed into three power spectrum densities called white, pink, and red (also know as Brownian). The noise density controls the amount of noise power per unit of bandwidth and it is proportional to 1/ f 0 (white), 1/ f (pink) and 1/ f 2 (red). This means that white 4

FSICR IA Super class is one class higher than DNV GL Ice(1A) according to Ref. [62], and design ice block thickness Hice is equal to 1.75 m. 14

Marine Structures 67 (2019) 102614

D. Polić, et al.

Fig. 14. Variation of the non-dimensional noise vector q (t ) due to the power spectrum density.

noise is equally powerful at any frequency, and pink and red noise have more power at the low frequencies (i. e pink noise have equal power in the frequency range from 1 to 10 Hz as in the band from 100 Hz to 1 kHz). Fig. 14 illustrates noise in all three power spectrum densities. The influence of white noise on the simulated propeller shaft torque, calculated from the PM model output n (t ) and using Eq. (1), is illustrated in Fig. 4. The same principle is applied for producing the artificial measured angular velocity ( ˆs ). 3.3. Model validation The dynamic response of the PM system exposed to torque Qp is calculated using the PM model and the Backward Euler (BE) solver with a fixed time step of 0.25 ms. Based on the torque Qp input to the propeller sub-model and combustion pressure pcyl input to the engine sub-model, the propeller shaft modal response n (t ) is calculated in the 20-sim-software version 4.4. Knowing the simulated response n (t ) , the location of the measurements xk and mode shape of the propeller shaft (Eq. (23)), one can obtain a time history of angular deformation s (Eq. (13)) and velocity s (the first time derivative of Eq. (3)). The angular deformation s (xk , t ) and velocity s (xk , t ) are then colored with the noise, set as 1% (e = 0.01)5 of the mean angular deformation ¯s = 2.38*10 4 rad (equivalent to 7.7 kNm) and velocity s = 9.34 rad/s, respectively. The obtained artificial measured angular deformation ˆs and velocity ˆs are then sampled with the four sampling frequencies that are listed in Table 4. In addition to different sampling frequencies, five different locations xk along the propeller shaft, given in Table 5, are considered as possible measurement locations. Together with four measurement setups, several possible combinations of inputs for each inverse approach model are created. The number of possible combinations per approach is given in Table 6. The cross-section of the propeller shaft is set by diameter Dshaft = 0.6 m and dshaft = 0.2 m. The colored and sampled deformation s (xk , t ) and angular velocity s (xk , t ) are then filtered using locally weighted scatter plot smoothing (loess MATLAB function). The loess function uses weighted linear least squares and a 2nd degree polynomial interpolation to reduce the random noise in the signal while a sharp step response is retained, see Fig. 4. The smoothing process is considered to be local because the ith data point that should be smoothed (called smooth value) is determined by neighbouring data points defined within a span. The span6 is optimized based on the inverse rule-based torque load. The optimal span for different sampling frequencies is given in Table 7. Furthermore, due to the requirement set by a system of liner equations, given in the inverse model, to have all input signals sampled with the same sampling frequency, signals collected with a frequency lower than 1 kHz are interpolated using the shapepreserving piecewise cubic interpolation (MATLAB syntax pchip7) to a time step of 1 ms. After coloring, sampling, filtering, and interpolating, the deformation s (xk , t ) and angular velocity s (xk , t ) are taken as inputs to the inverse model, and transformed into 2720 inverse propeller torque loads Qp* (t ) using five different model approaches. The number of modes considered in the inverse propeller shaft sub-model and in the PM model is the same 54, and hence, inaccuracy of the response calculated by the PM using 53 + modes is incorporated into the solutions. The two approaches a0 and a1 that do not consider the inverse propeller shaft sub-model, the response is corrected accordingly to Eq. (27). The performance of each inverse model approach is evaluated in terms of accuracy using residual Qp (t )

Qp (t ) = Qp* (t )

(29)

Qp (t ),

5 ± 1% of change in the angular velocity over the time step of 0.25 ms can cause change in the Qp* equal to 34 MNm, while the change of ± 1% in the angulr deformation can cause change in the Qp* equal to 0.6 MNm. 6

The span is constant over the whole smoothing process and it consists out of even or odd number of data points (the minimum number of data points is 5). At the beginning and at the end of the smoothing process when the number of neighbouring data points left or right from the smooth value is lower than the half of span, an asymmetric weight function is used. Otherwise, the weight function is symmetric. All data points inside the span have weight and influence on the fit (interpolation), but the data point that should be smoothed has the largest weight and the most influence on the fit. 7 The pchip interpolation method has a continuity level 1, one less than spline, and hence it oscillates less freely between the sample points. 15

Marine Structures 67 (2019) 102614

D. Polić, et al.

Table 4 List of different sampling setups and sampling frequencies in Hz used for measurement of angular deformation s (xk , t ) and velocity s (xk , t ) . Sampling setup

f

f

A B C D

1000 100 1000 1000

1000 1000 100 10

Table 5 Locations xk along the propeller shaft where angular deformation xk Lshaft

7 12

=

s

and velocity

s

are collected.

8 12

9 12

11 12

10 12

Table 6 The number of possible combinations of measured signals for each measurement setup (Fig. 2) by considering five possible locations along the shaft. The possibility of placing two or three sensors at the same location is not considered. Measurement setup

1

2

3

4

Approach No. of combinations

a0 5

a1 and a2 20

a3 10

a4 30

Table 7 Size of the span, as a function of the sampling frequency, used during signal processing. The sampling frequencies are in Hz. f , f

Span size

1000 100 10

84 8 5

and rainflow-counting algorithm. Considering that the value of the residual Qp (t ) is a function of time, see Fig. 15b, the likelihood of the residual Qp (t ) being within certain limits (i. e. ± 0.1Q¯ hyd ) is derived by integrating the area under the probability density function (PDF). The PDF and cumulative distribution (CDF) are illustrated in Fig. 15c and d, respectively. ˜ Qp is a normalized residual equal to

˜ Qp =

Qp (t ) , Q¯ hyd

(30)

where Q¯ hyd is the mean value of the torque Qp before ice-propeller interaction (85% of MCR) and equal to 770 kNm, see Fig. 16. This direct comparison of Qp* torque with Qp torque over time is suitable for evaluation of the mean value and phase shift between two signals, but may not capture sharp and short differences caused by the ice-propeller interaction (if a few short ice impacts are missed in the long time history by the inverse model, the CDF value will not be greatly affected). Hence, oscillations in the Qp and Qp* torque caused by the ice-propeller interaction are counted using rainflow-counting algorithm and illustrated as a function of the Qmax torque, see Fig. 15e. The influence of the location xk for each approach (i. e. 30 combinations for approach a4) is illustrated in the CDF and rainflow-counting histogram as an error bar, see Fig. 17. The bullet marker in the error bar indicates a mean value, while the hight of the error bar indicates the minimal and maximal value. Considering the CDF, rainflow-counting algorithm, and variability of both due to the location xk , the acceptable inverse torque Qp* should fulfil the following: - The mean CDF value must be greater than 95% for normalized residual ˜ Qp equal to 10%. - The variation of the CDF value must be less than ± 5% for ˜ Qp = 10%. - The amplitude of cycles observed in Qp* torque must be equal or grater than those observed in the original Qp torque. 16

Marine Structures 67 (2019) 102614

D. Polić, et al.

Fig. 15. Overview of the performance evaluation of the inverse model approaches.

Fig. 16. Start and end of the inverse model during rule-based ice-propeller interaction. 17

Marine Structures 67 (2019) 102614

D. Polić, et al.

Fig. 17. The inverse rule-based torque load Qp* (t ) time history derived from response that is sampled with setup A and using five different approaches a.

18

Marine Structures 67 (2019) 102614

D. Polić, et al.

Fig. 17. (continued)

19

Marine Structures 67 (2019) 102614

D. Polić, et al.

Fig. 17. (continued)

20

Marine Structures 67 (2019) 102614

D. Polić, et al.

Fig. 17. (continued)

21

Marine Structures 67 (2019) 102614

D. Polić, et al.

Fig. 18. The inverse rule-based torque load Qp* time history derived from response that is sampled with setup B and using approach a2.

22

Marine Structures 67 (2019) 102614

D. Polić, et al.

Fig. 19. The inverse rule-based torque load Qp* time history derived from response that is sampled with setup B and using approach a4.

23

Marine Structures 67 (2019) 102614

D. Polić, et al.

Fig. 20. The inverse rule-based torque load Qp* time history derived from response that is sampled with setup C and using approach a1.

24

Marine Structures 67 (2019) 102614

D. Polić, et al.

Table 8 List of approaches that managed to calculate acceptable Qp* load during rule-based ice-propeller interaction. Approaches, which satisfied the criteria related with the cycle amplitude, but failed to reach CDF value of 95%, are placed inside parenthesis (). Sampling setup

Noise without

A B C D

a2, a2, a3, a3,

a3, a4 a3, a4 a4 a4

white

pink

a2, a4 a2

red (a4) a2, (a4) (a4)

4. Results 4.1. Inverse rule-based torque load The results of transformation of propeller shaft angular deformation s and velocity s , obtained during the rule-based ice-propeller interaction, are illustrated in Figs. 17–20. The torque Qp* (t ) is calculated over a period of 6 s, starting at 280 s, approximately one second before the ice-propeller interaction. The endpoint is set by regaining back a steady angular velocity p (t ) . Fig. 16 illustrates the start and end points of the inverse model. As expected approach a0 fails to produce acceptable torque Qp* regardless of the presence of noise in the measured data, see Fig. 17 (the propeller load and propeller shaft response will never be similar in the transient state). The amplitudes of the half cycles in approach a0 are significantly smaller than expected and the CDF value is just 65% at ˜ Qp = 10%. Such derived inverse torque Qp* differs significantly from Qp , and therefore the performance of approach a0 is not discussed or illustrated further in the paper. Approaches a1-a4 give the acceptable inverse torque Qp* when s and s are sampled with setup A and are not subjected to the noise, see Fig. 17a. The performance of approaches a1-a4 in the presence of noise is illustrated in Fig. 17b-d. The presence of high-frequency noise (white and pink) in the measured data results in a smaller CDF value for smaller ˜ Qp , see Fig. 17b and c. Calculation of the inverse torque Qp* (t ) load using only two records of s (xk , t ) is unachievable. The torque difference between the two shaft ends, derived using a3 from two independent records of s (xk , t ) , is completely distorted in the presence of high-frequency noise and it can led to the error in the propeller angular acceleration s,2 (t )/ t . A small error in s,2 (t )/ t multiplied with the propeller inertia Jp leads to large under- or over-estimation of the Qp (t ) torque. The additional record of (xk , t ) , which is used to control and correct the sensitive s,2 (t )/ t , leads a4 to perform much better than a3. The CDF value and cycle amplitudes are more reasonable than those derived for a3, but still not good enough to be considered acceptable. Low-frequency noise (red) in the measured data, illustrated in Fig. 17d, is also very interesting because the red noise affects the Qp* load derived from a single record of * * s (xk , t ) differently than the Qp load derived from two independent s (xk , t ) records. The Qp cycle amplitudes derived using a1 and a2 are smaller than the Qp cycle amplitudes, whereas a4 has correct cycle amplitudes, but the mean value of the Qp* varies for each combination of the input signals, and hence, variation in the CDF value is significant in a4. The amplitude of the ice-propeller impact is more critical than the mean value for design, so the performance of a4 can be considered to be partially acceptable. Reducing the sampling frequency of angular deformation ( f ) from 1 kHz (setup A) to 100 Hz (setup B), the performance of a2 is moderately improved and the performance of a4 is worsened, see Figs. 18 and 19, respectively. In the presence of high-frequency

Fig. 21. Start and end of the inverse model during random ice-propeller interaction. 25

Marine Structures 67 (2019) 102614

D. Polić, et al.

noise (white and pink), the variation of the CDF and cycle amplitude is reduced in a2, whereas in a4 the latter is increased. Furthermore, the cycle amplitudes calculated by a2 in the presence of red noise are correct. The performance of other approaches is not affected by the change in the f frequency. Reducing the sampling frequency of angular velocity ( f ) from 1 kHz (setup A) to 100 Hz (setup C) does not have any positive effect on the performance of any approach. In the presence of high-frequency noise, a1 overestimates the amplitude and number of cycles in the Qp* load, see Fig. 20. Approaches a2 and a4 are unusable, and their performance is similar to a3 illustrated in Fig. 17b-c. Reducing f frequency further (setup D), the Qp* load calculated using approaches a2-a4 diverges from the Qp load, whereas in the case of a1, the Qp* load starts to be similar to Qp* derived using a0.

Fig. 22. CDF and number of half cycles observed in the inverse random torque Qp* (t ) , derived from response that is sampled with setup A and using four different approaches a. Approach a0 is omitted due to low accuracy in calculating the rule-based Qp* (t ) torque. 26

Marine Structures 67 (2019) 102614

D. Polić, et al.

Fig. 22. (continued)

The performance of all five inverse approaches, in the case of rule-based ice-propeller interaction, is summarized and presented in Table 8. The successful calculation of the acceptable Qp* load using a2 depends on the f frequency. Approach a3 and a4 are the only two approaches whose ability to calculate the acceptable Qp* load from a perfect measured signal does not depend on the sampling setup in the examined conditions. In the presence of noise in the measured data, limited success, calculating acceptable Qp* load, is achieved using a2 and a4. 4.2. Inverse random torque load Using the same rule for the start and end points as in the previous section, the torque Qp* during random ice-propeller interaction is calculated over a period of 1600 s. The start and end points of the inverse model are illustrated in Fig. 21, as well as the time history 27

Marine Structures 67 (2019) 102614

D. Polić, et al.

Fig. 23. Influence of discontinuities and filter (smoothing function) on the inverse random Qp* (t ) torque, derived from response that is sampled with setup A.

of the Qp torque and p angular velocity during random ice-propeller interaction. The random triangular definition of the ice impact, as opposed to the rule-based half-sinus definition with a fixed Cq and i ratio, results in lower performance of a1-a4 for all examined sampling setups. Fig. 22 illustrates the performance using sampling setup A. The cycle amplitude in the Qp* torque is almost always smaller than the corresponding one in the Qp torque for three reasons: a greater number of discontinuities in the random Qp torque time history, the presence of short and large ice impacts, and the size of the span 28

Marine Structures 67 (2019) 102614

D. Polić, et al.

Table 9 List of approaches that managed to calculate acceptable Qp* load during random ice-propeller interaction. Sampling setup

Noise without

A B C D

a1, a3, a3, a3,

white

pink

red

a3, a4 a4 a4 a4

considered during signal processing (filtering). Fig. 23 illustrates the influence of discontinuity and size of the span throughout the short segment of the inverse random torque Qp* time history with a single ice impact set by Cq 0.8 and i 15 . A spurious oscillation is calculated by the inverse approaches around the discontinuity point, and in the case of a2, the spurious oscillation exceeds the oscillation caused by the genuine ice impact, see Fig. 23a. The power of the spurious oscillation grows with the reduction of i size and with the reduction in the sampling frequency. By filtering the colored measured signals, the influence of discontinuity can be reduced as illustrated in Fig. 23b, but at the expense of accuracy because the sharp ice impacts are smoothed together with the discontinuity points. The results illustrated in Fig. 23 are derived by considering sampling setup A. The successful calculation of the acceptable Qp* , considering random ice impacts, for different sampling frequencies and noise colors is presented in Table 9. None of the approaches can successfully calculate Qp* when noise is present. Reducing the size of the span, compared to the value optimized for the rule-based ice-propeller interaction, helps to some extent. The difference in the Qp* and Qp cycle amplitudes during short and large ice impacts can be reduced, but the influence of the noise in the Qp* time history is greater in the case of high-frequency noise. The spurious oscillations calculated by the inverse approaches are frequent and independent from the discontinuity points and ice impacts. The low-frequency noise does not generate spurious oscillations, but does cause changes in the mean value. 4.2.1. Influence of starting time of the inverse model Considering two different starting points of the inverse model, before and during ice-propeller interaction, similar Qp* torque time histories are obtained, as illustrated in Fig. 24. The starting point of the inverse model during ice-propeller interaction is selected between two similar ice impacts, set by Cq 0.75 and i 85 . The only minor differences in the Qp* are observed in a4, and they are related with the correction of the calculated angular velocity s,2 (t ) using measured s (xk , t ). 5. Summary and discussion The procedure for calculating the propeller torque load based on the propeller shaft response is presented. In the absence of the model-scale or full-scale measurements, the used propeller shaft torque Qs (xk , t ) and angular velocity s (xk , t ) response are calculated using PM model, and colored with three types of artificial noise (white, pink and red) that introduced uncertainty. The proposed five inverse model approaches a0-a4, derived based on the four proposed measuring setups and three forms of system discretization, consider only torsional mechanisms of deformation in the elastic domain. Other mechanisms of deformation (e. g. axial and lateral) are not considered and therefore the presented procedure is more suitable for longer propeller shaft, traditionally used in the PM systems where a fixed or controlled pitch propeller is directly driven by the diesel engine with or without gearbox. The equation of torsional vibration is derived in the modal form by assuming angular position φ as a sum of products of timedepended generalized modal coordinates η and kinematically valid mode shapes Y. Consequently, each mode is described with an own equation of motion. In the inverse model, the equation of motion of flexible modes is simplified by neglecting modal inertia and modal damping term. Beside the propeller shaft and a rigid propeller, the rest of the PM system is neglected in the inverse model. The inverse propeller torque Qp* (t ) load is derived by inverse approaches and compared with the corresponding propeller torque Qp (t ) load (input to the PM model). Considering the residual (a difference between Qp* and Qp torque) and counting the oscillation cycles in both loads, it is possible to conclude: 1. Propeller shaft torque response Qs (xk , t ) cannot be considered to be equal or similar to the propeller torque Qp (t ) during the icepropeller interaction - derived from a0 2. In the presence of noise in the propeller shaft response data, limited success, calculating Qp* (t ) torque, is achieved. The accuracy of torque Qp* (t ) is generally greater if propeller shaft angular velocity response s (xk , t ) is utilised in the inverse approach - derived from a1, a2 and a4 3. Calculation of the propeller angular acceleration ¨s,2 (t ) , and hence torque Qp* (t ) , from the two independent propeller shaft torque responses Qs (xk , t ) is vary susceptible to the presence of noise - derived from a3 4. If the number of measured signals is greater than the number of inputs required by the inverse propeller shaft sub-model, additional signal(s) can be used to increase the accuracy of the Qp* (t ) torque by controlling and correcting the output from the inverse propeller shaft sub-model - derived from a4 5. sampling frequency f is equally important as f frequency for calculation of Qp* (t ) torque - derived from a1, a2 and a4 6. It is easier to calculate the inverse propeller torque Qp* (t ) load during the rule-based short sequence of half-sinus ice-propeller 29

Marine Structures 67 (2019) 102614

D. Polić, et al.

Fig. 24. Influence of starting time of the inverse model, in the presence of white noise, on the inverse random Qp* (t ) torque. Torque Qp* (t ) is derived from response that is sampled with setup A and colored with pink noise.

impacts than Qp* (t ) load during the long sequence of occasional and random triangular impacts 7. A segment of the Qp* (t ) time history with short and large ice impacts proves to be the most challenging to calculate because it is easy to partially or completely filter out such impact from the measured signal 8. The calculated Qp* (t ) time history is not dependent on starting time of the inverse approach 30

Marine Structures 67 (2019) 102614

D. Polić, et al.

9. All inverse approaches gave well-posed Qp* (t ) torque solution regardless of the shape of Qp* (t ) torque, presence of noise in the data, measurement setup and location, sampling frequency, and starting time of the inverse approach 10. The Qp* (t ) torque can be calculated at a certain time instant based on the two points in the measured propeller shaft response sequence. The derived conclusions are based on the assumption that noise amplitude is proportional to 1% of mean angular deformation and velocity. Because 1% of the change in the angular deformation and velocity over the same time step does not have the same impact on the Qp* load (noise in the angular velocity has approximately 55 times more influence that noise in the angular deformation), the conclusions may be bias towards approach a4. Approach a4 considers two independent propeller shaft deformations s (xk , t ) as main input to inverse propeller shaft sub-model and uses propeller shaft angular velocity s (xk , t ) to control the output from the inverse propeller shaft sub-model. Further development of control used in the approach a4 would be needed. The second most promising approach, beside a4, is a2, in the presence of noise. Approach a2 considers propeller shaft deformations s (xk , t ) and angular velocity s (x k , t ) as main input to inverse propeller shaft sub-model and does not have any control over the output. The artificial noise added to the calculated propeller shaft response may differ from the real noise, and hence, the proposed procedure must be validated with real data and real noise that are obtained from a model-scale test or full-scale trial. The model-scale test or full-scale trial must allow a simultaneous measurements of the propeller torque load and propeller shaft response in order to provide sufficient amount of data for validation. After validation, a successful use of inverse methods would enable shaft line measurements to be used in the development of the ice-propeller prediction models and classification regulations concerning design torque for the PM system. Sensitivity of the each approach to measured signal location along the propeller shaft is illustrated as a part of the approach performance evaluation, but without evaluation of the influence of permutation of the five location and ten noise vectors on the Qp* load, it is not possible to suggest the best location to measure angular deformation or velocity. Acknowledgment This research is funded through the Norwegian Research Council Project no. 194529. References [1] Carlton JS. Propulsion systems. third ed. Oxford: Butterworth-Heinemann; 2012. [2] J. W. Lewis, F. W. DeBord, V. A. Bulat, Resistance and propulsion of ice-worthy ships, Society of Naval Architects and Marine Engineers 90 (8). [3] Jussila M. Ice loads on the propulsion system of an ice breaking tug. 7th international conference on port and ocean engineering under arctic conditions, vol. 2, helsinki, Finland. 1983. p. 575–90. [4] Laskow V, Revill C. Engineering analysis of ice/propeller interaction data, Tech. Rep. TP 8450E. Transportation Development Centre Report; 1986. [5] Kannari P. Full scale and model tests performed with a nozzle and open propeller simultaneously. 9th IAHR internation symposium on ice. Japan: Sapporo; 1988. p. 772–81. [6] Jussila M, Koskinen P. Ice loads on propeller blade of small car ferry. 10th international conference on port and ocean engineering under arctic conditions. 1989. p. 862–72. [7] Jussila M, Koskinen P. Ice loads on cp-propeller and propeller shaft of small ferry and their statistical distributions during winter ’87. 8th international conference on offshore mechanics and arctic engineering, vol. 4. Hague: Netherland; 1989. p. 351–9. [8] Keinonen AJ, Browne R. Ice propeller interaction forces, Tech. Rep. TP10401E. Transportation Development Centre Report; 1990. [9] Cowper B, Browne R, Glen I, Ritch R. Resistance and propulsive performance trials of the mv terry fox and mv ikaluk in level ice. Trans - Soc Nav Archit Mar Eng 1992;100:315–43. [10] Williams F, Spencer D, Mathews S, Bayly I. Full scale trials in level ice with canadian r-class icebreaker. Trans - Soc Nav Archit Mar Eng 1992;100:293–313. [11] Wang J. Prediction of propeller performance on a model podded propulsor in ice (propeller-ice interaction) PhD. thesis St. John’s, Canada: Faculty of Engineering and Applied Science, Memorial University of Newfoundland; April 2007. [12] Searle S, Veitch B, Bose N. Experimental investigation of a highly skewed propeller in ice. J Offshore Mech Arct Eng 2001;123(4):191–7. [13] Moores C, Veitch B, Bose N, Jones S, Carlton J. Multi-component blade load measurements on a propeller in ice. Trans - Soc Nav Archit Mar Eng 2002;110:169–87. [14] Wang J, Akinturk A, Jones SJ, Bose N, Kim M-C, Chun H-H. Ice loads acting on a model podded propeller blade. J Offshore Mech Arct Eng 2007;129(3):236–44. [15] Karulina M, Karulin E, Belyashov V, Belov I. Assessment of periodical ice loads acting on screw propeller during its interaction with ice. International conference and exhibition on ships and structures in ice, vol. 163. 2008. p. 1–7. [16] Hagesteijn G, Brouwer J, Bosman R. Development of a six-component blade load measurement test set-up for propeller-ice impact. International conference on ocean, offshore and arctic engineering, vol. 84192. 2012. p. 1–10. [17] von Bock Und RUF, Ehlers S, et al. On the scalability of model-scale ice experiments. J Offshore Mech Arct Eng 2015;137(5):051502. [18] von Bock Und RUF. The mechanical behavior of model-scale ice: experiments, numerical modelling and scalability PhD. thesis Helsinki, Finland: School of Engineering, Aalto University; June 2016. [19] Browne RP, Revill CR, Keinonen AJ. Propeller design load model, Tech. Rep. TP 13243E. Institute for Ocean Technology, National Research Council of Canada; April 1998. [20] Ikonen T, Peltokorpi O, Karhunen J. Inverse ice-induced moment determination on the propeller of an ice-going vessel. Cold Reg Sci Technol 2015;112:1–13. [21] Wind J. The dimensioning of high power propeller system for arctic icebreaking vessels. 5th lips propeller symposium. Drunen; 1983. p. 1–30. [22] Kotras T, Humphreys D, Baird A, Morris G, Morley G. Determination of propeller ice milling loads. J Offshore Mech Arct Eng 1985;2:336–43. [23] Chernuka M, Jategaonkar R, Norwood M, Warner J. Development of a procedure for predicting propeller-ice interaction forces, Tech. Rep. TP 9850E. Transport Canada Publication; 1989. [24] Veitch B. Predictions of ice contact forces on a marine screw propeller during the propeller-ice cutting process PhD. thesis Helsinki, Finland: Acta Polytechnica Scandinavica; 1995. [25] Soininen H. A propeller-ice contact model PhD. thesis Espoo, Finland: Helsinki University of Technology; 1998. [26] Ye L, Wang C, Chang X, Zhang H. Propeller-ice contact modeling with peridynamics. Ocean Eng 2017;139:54–64. [27] Koskikivi J, Kujala P. Long-term measurements of ice induced loads on the propulsion machinery of product tanker sotks, Tech. Rep. vol. 42. Espoo, Finland: Technical Research Centre of Finland (VTT); October 1985. [28] Browne RP. Analysis of canadian full scale propeller and ice interaction trials data for correlation with empirical models. Tech. Rep. CR-1997-12. Institute for

31

Marine Structures 67 (2019) 102614

D. Polić, et al.

Ocean Technology, National Research Council of Canada; 1997. [29] Dahler G, Stubbs JT, Norhamo L. Propulsion in ice – big ships. 9th international conference and exhibition on performance of ships and structures in ice, vol. 133. 2010. p. 1–10. [30] Bekker A, Suominen M, Peltokorpi O, Kulovesi J, Kujala P, Karhunen J. Full-scale measurements on a polar supply and research vessel during maneouver tests in an ice field in the baltic sea. International conference on ocean, offshore and arctic engineering, no. 24128. American Society of Mechanical Engineers; 2014. p. 1–10. [31] de Waal R, Soal K, Bekker A, Heyns P. Rotational dynamic characteristics of the shaft line of a rotational dynamic characteristics of the shaft line of a polar supply and research vessel during open water, cavitation and ice navigation. ISMA 2016 - international conference on noise and vibration engineering. 2016. p. 1–15. [32] de Waal R, Bekker A, Heyns P. Bi-polar full-scale measurements of operational loading on polar vessel shaft-lines. POAC 2017 - international conference on port and ocean engineering under arctic conditions, no. 74. 2017. p. 1–18. [33] Jones SJ, Soininen H, Jussila M, Koskinen P, Newbury S, Browne R. Propeller-ice interaction. Trans - Soc Nav Archit Mar Eng 1997;105:399–425. [34] He Q, Du D. Modeling and calculation analysis of torsional vibration for turbine generator shafts. J Inf Comput Sci 2010;7:2174–82. [35] Sarkar TK, Tseng FI, Rao SM, Dianat SA, Hollmann BZ. Deconvolution of impulse response from time- limited input and output: theory and experiment. IEEE Transactions on Instrumentation and Measurement 1985;34(4–1):541–6. [36] Jacquelin E, Bennani A, Hamelin P. Force reconstruction: analysis and regularization of a deconvolution problem. J Sound Vib 2003;265(1):81–107. [37] Sanchez J, Benaroya H. Review of force reconstruction techniques. J Sound Vib 2014;333(14):2999–3018. [38] Polić D, Æsøy V, Ehlers S. Transient simulation of the propulsion machinery system operating in ice - modeling approach. Ocean Eng 2016;124:437–49. [39] Cleveland WS. Robust locally weighted regression and smoothing scatterplots. J Am Stat Assoc 1979;74(368):829–36. [40] Polić D, Ehlers S, Æsøy V. Propeller torque load and propeller shaft torque response correlation during ice-propeller interaction. J Mar Sci Appl 2017;16(1):1–9. [41] Polić D, Ehlers S, Æsøy V, Pedersen E. Shaft response as a propulsion machinery design load. 33rd international conference on ocean, offshore and arctic engineering, no. 23880, san francisco, USA. 2014. p. 1–10. [42] Margolis DL. A survey of bond graph modelling for interacting lumped and distributed systems. J Frankl Inst 1985;319(1):125–35. [43] Meirovitch L. Analytical methods in vibrations vol. 438. New York: Macmillan; 1967. [44] Daniel JI. Engineering vibration. third ed. New Jersey: Prentice-Hall, Inc.; 2001. [45] Rao SS. Vibration of continuous systems. John Wiley & Sons; 2007. [46] Karnopp DC, Margolis DL, Rosenberg RC. System dynamics - modeling and simulation of mechatronic systems. fourth ed. John Wiley & Sons, Inc.; 2005. [47] Borutzky W. Bond graph methodology: development and analysis of multidisciplinary dynamic system models. first ed. Incorporated: Springer Publishing Company; 2010. [48] Borutzky W, editor. Bond graph modelling of engineering systems - theory, applications and software support. first ed.NY, NY, U.S.A: Springer-Verlag; 2011. [49] Karnopp DC. Computer representation of continuous vibratory systems using normal modes and bond graph techniques, in: soc. Comput. Simul. 1968. p. 129–35. [50] Karnopp DC, Rosenberg RC. Analysis and simulation of multiport systems - the bond graph approach to physical system dynamics. Journal of the Franklin Institute 1968;286(6). [51] Engja H, Bunes O. Modelling and simulation of propulsion dynamics due to clutching. Proceedings of the 6th international symposiumon marine engineering, no. 25, tokyo,Japan. 2000. p. 151–8. [52] Bruun K, Pedersen E, Valland H. Modelling for transient torsional vibration analysis in marine power train systems. Proceedings of the 2005 international conference on bond graph modeling. 2005. p. 185–90. [53] Polić D, Ehlers S, Æsøy V, Pedersen E. Propulsion machinery operating in ice - a modelling and simulation approach. European Council for modelling and simulation, alesund, Norway. 2013. p. 1–7. [54] Taskar B, Yum KK, Steen S, Pedersen E. The effect of waves on engine-propeller dynamics and propulsion performance of ships. Ocean Eng 2016;122(Supplement C):262–77. [55] Yum KK, Taskar B, Pedersen E, Steen S. Simulation of a two-stroke diesel engine for propulsion in waves. International Journal of Naval Architecture and Ocean Engineering 2017;9(4):351–72. [56] Wilson WK. Practical solution of torsional vibration problems : 1. Frequency calculations. third ed. London: Chapman & Hall; 1963. [57] Batrak Y, Serdjuchenko A, A.I. T. Calculation of propulsion shafting transient torsional vibration induced by ice impacts on the propeller blades. Proceedings of world maritime technology conference, Russia, saint-petersburg. 2012. p. 1–8. [58] Batrak Y, Serdjuchenko A, T.A.I.. Calculation of torsional vibration responses in propulsion shafting system caused by ice impacts. Proceedings of 1st torsional vibration symposium. Austria: Salzburg.; 2014. p. 1–8. [59] Persson S. Ice impact simulation for propulsion machinery. Proceedings of the 1st torsional vibration symposium. Austria: Salzburg; 2014. p. 34–41. [60] DNV GL (Det Norske Veritas, Germanischer Lloyd. DNVGL-RU-SHIP Pt.6 Ch.6 Cold climate. July 2018. [61] IACS (International Association of Classification Societies). UR I3 Machinery requirements for polar class ships. 2007. [62] FSICR (Finnish Swedish Ice Class Rules). Ice class regulations. 2010. [63] Veitch B. Propeller-ice interaction PhD. thesis Espoo, Finland: Helsinki University of Technology; 1992. [64] Wirgin A. The inverse crime. ArXiv Mathematical Physics e-prints; 2004. [65] Kaipio J, Somersalo E. Statistical and computational inverse problems vol. 160. Springer Science & Business Media; 2006.

32