2.1 Dynamic analysis of diesel engine piston

Taking the force and motion of a single cylinder of a diesel engine as an example, as shown in Figure 1, under the condition of healthy and uniform speed, the single-cylinder crank rotates at equal angular speed, the crank length is R, the connecting rod length is L, and the crank moves to the highest point and the lowest point corresponding to the top dead centre and bottom dead centre of the cylinder piston, respectively. Therefore, the dynamic equation is established based on the displacement of piston motion as

(1)\begin{equation}x = A{A_1} = R + L - L\cos \alpha - R\cos \beta\end{equation}

Figure 1.

Piston model

According to the triangle △OAB

(2)\begin{equation}L\textrm{sin}\alpha \textrm{ = }R\textrm{sin}\beta\end{equation}

Defining $\lambda \textrm{ = (}R/L\textrm{) = (sin}\alpha /\textrm{sin}\beta \textrm{)}$, the above formula can be expanded by Newton binomial to obtain the following form:

(3)\begin{equation}\textrm{cos}\alpha \textrm{ = 1} - \textrm{0}\textrm{.5}{\lambda ^\textrm{2}}\textrm{si}{\textrm{n}^\textrm{2}}\beta\end{equation}

By substituting the original formula, the piston displacement formula is

(4)\begin{align} x & =A{A_1} = R + L - L(1 - 0.5{\lambda ^2}{\sin ^2}\beta ) - R\cos \beta \\ & =R(1 - \cos \beta ) + 0.5R\lambda {\sin ^2}\beta \end{align}

The v (velocity) and a (acceleration) of the x (piston) can be obtained by calculating the first derivative and the second derivative of the displacement, respectively:

(5)\begin{equation}\left\{\begin{array}{@{}l} v = \dot{x} = R\omega (\sin \beta + 0.5\lambda \sin 2\beta )\\ a = \ddot{x} = R{\omega^2}(\cos \beta + \lambda \cos 2\beta ) \end{array} \right.\end{equation}

The acceleration vibration signal of the single-cylinder piston of a diesel engine is essentially an angle cyclostationary signal. During operation, the inertia force generated by the piston during periodic reciprocating motion and the gas pressure in the cylinder are the main excitation sources of diesel engine cylinder vibration. The inertial force can be expressed as

(6)\begin{equation}{P_j} = {m_j}a = {m_j}R{\omega ^2}(\cos \beta + \lambda \cos 2\beta ) = {m_j}R{\omega ^2}\cos \beta + {m_j}R{\omega ^2}\lambda \cos 2\beta\end{equation}

where ${m_j}$ is the mass of the moving part of each cylinder, and the primary inertia force and secondary inertia force can be expressed as

(7)\begin{equation}\left\{\begin{array}{@{}l} {P_j}_1 = {m_j}R{\omega^2}\cos \beta \\ {P_j}_2 = {m_j}R{\omega^2}\lambda \cos 2\beta \end{array} \right.\end{equation}

For an odd-cylinder diesel engine, the first-stage inertia force will be partially offset. In contrast, the second-stage inertia force will be superimposed, and the amplitude-doubling effect will be displayed in the spectrum. The cylinder pressure of a diesel engine can be expressed as

(8)\begin{equation}P = {P_g}\frac{{\pi {D^2}}}{4}\end{equation}

where ${P_g}$ is the cylinder pressure and D is the cylinder diameter. According to the operating characteristics of the four-stroke diesel engine, the crankshaft produces a pressure peak every two cycles, so it will produce 0 ⋅ 5 times the crankshaft rotation frequency.

The total piston force can be expressed as

(9)\begin{equation}\sum {P = P + {P_j} = {P_g}\frac{{\pi {D^2}}}{4} + {m_j}R{\omega ^2}\cos \beta + {m_j}R{\omega ^2}\lambda \cos 2\beta }\end{equation}

The torsional vibration force transmitted to the shafting can be obtained by synthesis and decomposition:

(10)\begin{equation}M = \sum P \ast \tan \alpha \ast (L\cos \alpha + R\cos \beta )\end{equation}

In practical engineering, the structure and motion of a diesel engine lead to a complex and changeable dynamic system which will not only be affected by these internal forces but also external factors. Therefore, the measured diesel engine vibration is a multi-component complex non-stationary signal.

Remark 1: Due to the complex structure and diversity of the marine environment, evaluation standards have also become problematic in scientific research and engineering applications. The evaluation of the vibration of mechanical equipment or systems has always been a hot topic of study. The vibration evaluation of a simple mechanism might be well known, but the influencing factors of vibration of large and complex equipment are not mature or well understood. Due to the complexity of the marine environment and many kinds of vibration sources and noise during actual operations, the vibration signal of a Ro-Ro passenger ship is highly complex. Extracting the single component with the most extensive feature and removing the interference signal from the complex signal is a difficult problem in engineering.

Remark 2: The different multi-model vibration signals generated by the Ro-Ro passenger ship's sensors are a significant challenge for practical engineering analysis.

2.2 Iterative generation differential decomposition

The generated differential equation can be regarded as a mapping transformation of a signal (Li and Feng, Reference Li and Feng2017). This transformation can decompose the signal into a set of solutions of a series of differential equations to estimate the signal energy, amplitude, frequency and other information. Without losing generality, a harmonic oscillation signal can be expressed in the following form:

(11)\begin{equation}x(t) = A(t)\cos (\omega t + \varphi )\end{equation}

where $A(t)$ is the amplitude modulated signal, $\omega$ is the natural angular frequency and $\varphi$ is an initial phase.

First, define the following basic state functions of $x(t)$:$\delta x$ is the relative rate of change representing the function $\dot{x}(t)\textrm{ = } - A(t)\omega \sin (\omega t + \varphi )$, $\delta \textrm{y}$ is the relative rate of change of the representative function $\ddot{x}\textrm{ = } - A(t){\omega ^2}\cos (\omega t + \varphi )$, and ${k_x}$ and ${\lambda _x}$ are the derivative of $\delta x$ and $\delta \textrm{y}$:

(12)\begin{equation}\left\{ \begin{array}{@{}l} {\delta_x} = \dfrac{{\dot{x}}}{x} = \dfrac{{dx/x}}{{dt}}\\ {\delta_y} = \dfrac{{\dot{y}}}{y} = \dfrac{{\ddot{x}}}{{\dot{x}}} = \dfrac{{dy/y}}{{dt}}\\ {k_x} = {\delta_x}{\delta_y} = \dfrac{{\ddot{x}}}{x}\\ {\lambda_x} = \dfrac{{{\delta_y}}}{{{\delta_x}}} = \dfrac{{dy/y}}{{dx/x}} \end{array} \right.\end{equation}

The differential equation can be deduced from the first-order differential and second-order differential of the signal $x(t)$ as

(13)\begin{equation}\ddot{x} - [{k_A} - {\delta _A}(2{\delta _A} + {\delta _\omega }) - {\omega ^2}]x - (2{\delta _A} + {\delta _\omega })\dot{x} = 0\end{equation}

where ${\delta _A} = \dot{A}/A$, ${\delta _\omega } = \dot{\omega }/\omega$, ${k_A} = \ddot{A}/A$.

Generally, the modulation part $A(t)$ and $\varphi (t)$ changes slowly with respect to the carrier signal, so they can be approximately constant. Therefore, the state function of amplitude and frequency can be approximated ${\delta _A} = 0$, ${\delta _\omega } = 0$, ${K_A} = 0$. Then, the solvable instantaneous frequency and envelope amplitude are

(14)\begin{equation}\left\{ \begin{array}{@{}l} {\omega^2} ={-} \dfrac{{\ddot{x}}}{x} ={-} {k_x} ={-} {\delta_x}\delta (x{\delta_x})\\ {A^2} = {x^2} - \dfrac{{x{{\dot{x}}^2}}}{{\ddot{x}}} = {x^{_2}}\left( {1 + \dfrac{1}{{{\lambda_x}}}} \right) \end{array} \right.\end{equation}

It is necessary to replace the differential of continuous signals with a central difference and the derivative with a central difference quotient for analysing discrete signals:

(15)\begin{equation}\left\{ \begin{array}{@{}l} \dot{x}(n) = \dfrac{{x(n + 1) - x(n - 1)}}{{2\Delta t}}\\ \ddot{x}(n) = \dfrac{{x(n + 2) - 2x(n) + x(n - 2)}}{{4\Delta {t^2}}} \end{array} \right.\end{equation}

The estimated values of the instantaneous frequency ω and envelope amplitude A of the discrete signal can be obtained through Equations (14) and (15). By using this algorithm, not only can the envelope be tracked, but also any drift in the frequency can be tracked simultaneously

(16)\begin{equation}\left\{ \begin{array}{@{}l} \omega = \sqrt {\dfrac{{2x(n) - x(n + 2) - x(n - 2)}}{{4x(n)}}} \\ A = \sqrt {{x^2}(n) - \dfrac{{{{(x(n + 1) - x(n - 1))}^2}}}{{x(n + 2) - 2x(n) + x(n - 2)}}x(n)} \end{array} \right.\end{equation}

Without loss of generality, the multi-component vibration signal can be expressed as the superposition of multiple AM-FM signals as

(17)\begin{equation}x(t) = \sum\limits_{i = 0}^M {a_i}(t)\cos [\int{\omega _i}(t){\rm d}t]\end{equation}

where ${a_i}(t)$ and ${\omega _i}(t)$ are the amplitude envelope and carrier frequency of the i-th component.

Demodulate the generated differential equation of Equation (17):

(18)\begin{equation}\left\{ \begin{array}{@{}l} {\omega^2} ={-} {k_x} = {\omega_1}^2\textrm{ + }{R_{\omega (t)}}\\ {A^2} = {x^2} - \dfrac{{x{{\dot{x}}^2}}}{{\ddot{x}}} = {A_1}^2(t) + {R_{A(t)}} \end{array} \right.\end{equation}

where ${R_{\omega (t)}}$ and ${R_{A(t)}}$ are the frequency and amplitude of the residual component, respectively. Because the residual component has high-frequency and low-energy characteristics, so it can be further separated by filtering.

For some specific reasons, such as the multi-component signals to be decomposed are discrete signals with a finite length, the zero-phase filter will produce distortion when eliminating phase shift at the initial and final ends of the signal, resulting in an endpoint effect. In addition, with the continuous increase of iterative operation, this ‘distortion’ will continue accumulating and cause ‘distortion diffusion’. Moreover, because the iterative differential equation decomposition method uses both difference and difference quotients to process and analyse discrete signals, outliers may appear at the endpoints. Consequently, a signal continuation method can effectively solve the problem of the end effect. The extension signal $\tilde{x}$ is obtained by the following symmetrical extension:

(19)\begin{equation}\tilde{x} = [x({N_e} - {N_e}_i + 2)\quad x\quad x(N - {N_e}_i)]\end{equation}

where N is the length of the original signal and ${N_e}$ is the extension length ${N_{ei}} = 1,2, \cdots ,{N_e}$. After the continuation, the signal will be extended, so after the first filtering, it is necessary to cut the continuation signal to restore the original length, and carry out the continuation again before the next iteration until the end of the iteration. The process flow of our proposed method is shown in Figure 2.

Figure 2.

Procedure

The main highlight of the proposed strategy lies in that the iterative-generation differential decomposition overcomes of the drawbacks of multi-component vibration signals.

3.1 Test experiment

The vibration test object of ship engineering is a diesel engine group of a specific type of Ro-Ro passenger ship. The single group engine is a 4-stroke, non-reversible, turbocharged, inter-cooled diesel engine with direct fuel injection, which directly injects fuel to output power through nine cylinders, with a rated speed of 1200 rpm and a maximum continuous output rated power of 1665 kW.

According to the field test requirements, the measuring points are arranged in the side spiral tube (Point 1), the middle of the front chassis (Point 2) and the bottom-corner chassis (Point 3). The measuring directions are in the horizontal diameter direction, horizontal axial direction and vertical direction, respectively. The model of signal acquisition instrument is a Bruel & Kjaer 3505, and PCB acceleration sensors greatly improve the performance of the tests. Three working conditions with loads of 0%, 75% and 100% at rated speed are tested. There is no other interference in the field test environment, and the test site is shown in Figure 3.

Figure 3.

Test environment

3.2 Time-domain analysis

1. (1) Time domain analysis is used mainly to analyse the stability, transient and steady-state performance of the system by extracting time-domain characteristics. Therefore, firstly, we compare and analyse different load conditions at the same measuring point. The analysis results are shown in Figure 4. It can be seen that the dimensional eigenvalues of measuring Point 1, including kurtosis, skewness, variance, standard deviation and effective value, increase with an increase of load, which intuitively shows the energy change of signal. The dimensionless eigenvalues include waveform index, peak index, pulse index and margin index, which change irregularly with an increase of load. This is because the dimensionless index is highly sensitive to impact, but its stability is insufficient. In practical analysis, the dimensionless eigenvalues can be used to judge the energy distribution of a system operation. The dimensionless index is used to judge the system impact degree under the same working condition in different stages.

2. (2) The power spectrum shows that the power energy is concentrated on 360 Hz, 720 Hz, 1080 Hz  … corresponding to $n \times 4 \times 9 \times 0.5 \times {f_0}$, respectively, which is also in line with the vibration characteristics of the 4-stroke and 9-cylinder diesel engine of the test object. Moreover, with an increase of load, the amplitude of power excitation energy increases gradually.

3. (3) It can be seen that there is little difference in vibration amplitude under the three different loads, but in the corresponding power spectrum, the corresponding power excitation capacity of the diesel engine in the high-frequency region increases with an increase of load.

Figure 4.

Time domain analysis and comparison of three load conditions: (a) measuring Point 1; (b) measuring Point 2; (c) measuring Point 3

3.3 Frequency domain analysis

3.3.1 Iterative-generating differential decomposition

Firstly, the measured data are cleaned and sorted, and then the stable data of 10 s is intercepted for the next analysis. The first step is to generate differential decomposition for the data of three measuring points, as shown in Figure 5. To facilitate comparative analysis, the blue line represents the signal of measuring Point 1, the red line represents the signal of Point 2 and the green line represents the signal of Point 3. To reduce the proportion of the diagram, the first four decomposed components are selected for display.

Figure 5.

Iterative generating differential decomposition results

To better extract the signal dynamic parameters, selecting the best sensitive component for the next analysis is necessary. Firstly, the correlation coefficient between each component and the original signal is calculated, and the correlation threshold is set to $\lambda = \max ({\mathop{\rm cov}} f(i))/\beta$. It can be seen that the signal correlation of the three measuring points is the first component in Figure 6. Therefore, the first component is selected as the sensitivity analysis for the next analysis.

Figure 6.

The correlation coefficient

3.3.2 Fourier spectrum analysis

Fourier spectrum analysis is carried out for the first component decomposed by measuring Point 1. As shown in Figure 7, in the low-frequency region, 0 ⋅ 5-order rotation frequency and 3 times rotation frequency are not prominent, while 4 ⋅ 5-order rotation frequency is quite prominent. This is because the diesel engine is a 9-cylinder single-row cylinder, so the frequency conversion excitation 9 × 0 ⋅ 5 × ${f_0}$ is prominent. Other rotation frequencies and their harmonics, such as 3 ⋅ 5-order rotation frequency, 4 ⋅ 5-order rotation frequency, etc., are also prominent. It is usual for rotation frequency harmonics to occur in low-frequency areas. The relative distribution of healthy diesel engines is also relatively uniform. As shown in Figure 7, the overall amplitude increases with an increase in load, especially with the 4 ⋅ 5-order frequency.

Figure 7.

Fourier spectrum of measuring Point 1 under three load conditions: (a) Fourier spectrum; (b) comparison of results

Fourier spectrum analysis is carried out for the first component decomposed by measuring Point 2. As shown in Figure 8, compared with Point 1, the overall amplitude increases significantly in the low-frequency region. Point 2 is close to the cylinder head, and the vibration impact is more evident than in other positions. In terms of frequency conversion and frequency doubling, the prominent frequency amplitude is similar to that of Point 1, and the overall amplitude increases with the increase of load, especially the rotation frequency harmonics 4 ⋅ 5${f_0}$.

Figure 8.

Fourier spectrum of measuring Point 2 under three load conditions: (a) Fourier spectrum; (b) comparison of results

As shown in Figure 9, Fourier spectrum analysis is performed on the first component decomposed by measuring Point 3. The overall amplitude increases with an increase in load, where the performance of rotation frequency harmonics 1 ⋅ 5${f_0}$, 3${f_0}$ and 4 ⋅ 5${f_0}$ of 0% load is evident. In contrast, the amplitude of rotation frequency harmonics 1 ⋅ 5${f_0}$, 3${f_0}$ and 4 ⋅ 5${f_0}$ of 75% load increases significantly, and the amplitude of frequency at the same position of 100% load increases more significantly. It shows a significant incremental effect.

Figure 9.

Fourier spectrum of measuring Point 3 under three load conditions: (a) Fourier spectrum; (b) comparison of results

Table 1 details the statistical information of all Fourier spectrum analysis of three different loads at three measuring points, including frequency values and amplitudes of 0 ⋅ 5-order, 1 ⋅ 0-order, 1 ⋅ 5-order, 3 ⋅ 0-order and 4 ⋅ 5-order.

Table 1.

Fourier spectrum amplitude statistics

3.3.3 Envelope spectrum analysis

Firstly, envelope spectrum analysis is carried out for the first component decomposed by measuring Point 1. As shown in Figure 10, the amplitude of the envelope spectrum after demodulation significantly changes compared with that of the Fourier spectrum. The overall amplitude is higher, and the amplitudes of rotation frequency harmonics 1 ⋅ 5${f_0}$, 3${f_0}$ and 4 ⋅ 5${f_0}$ are higher. After demodulation, the resonance interference component is removed. Therefore, the bottom performance of the envelope spectrum is more average than that of the Fourier spectrum. The overall amplitude does not change significantly with load. However, the rotation frequency harmonics 4 ⋅ 5${f_0}$ increases significantly with an increase in load.

Figure 10.

Envelope spectrum of measuring Point 1 under three load conditions: (a) envelope spectrum; (b) comparison of results

Continuing to analyse measuring Point 2, as shown in the Figure 11, the amplitudes of 1 ⋅ 5-order rotation frequency, 3-order rotation frequency and 4 ⋅ 5-order rotation frequency are higher, the overall spectrum performance is more average and the characteristic frequency is more prominent than that of Point 1. And the overall amplitude changes obviously with the load. This is because Point 2 is close to the cylinder head, and strong vibration and a short transmission path make the analysis results more accurate and reliable.

Figure 11.

Envelope spectrum of measuring Point 2 under three load conditions: (a) envelope spectrum; (b) comparison of results

Envelope spectrum analysis is carried out for Point 3. As shown in Figure 12, the amplitude of 0 ⋅ 5-order rotation frequency, 1 ⋅ 5-order rotation frequency and 4 ⋅ 5-order rotation frequency is higher, and the overall amplitude of 3-order frequency is lower. This is because Point 3 is close to the foundation, and part of the energy of even times of frequency conversion impact is absorbed by the ship plate during transmission, resulting in lower amplitude. The overall trend is the same as that for other measuring points, which shows an increasing trend with the increase of load.

Figure 12.

Envelope spectrum of measuring Point 3 under three load conditions: (a) envelope spectrum; (b) comparison of results

Table 2 details the statistical information of all envelope spectrum analysis of three different loads at three measuring points, including frequency values and amplitudes of 0 ⋅ 5-order, 1 ⋅ 0-order,1 ⋅ 5-order, 3 ⋅ 0-order and 4 ⋅ 5-order.

Table 2.

Envelope spectrum amplitude statistics