Research on Vibration Excitation Estimation of Warship Propulsion Shafting Based on Waveguide Theory

The propeller excitation can be estimated by measuring the vibration response of the shafting, but there are many difficulties, such as difficult to implement, difficult to determine the boundary conditions and so on. For this reason, the waveguide model of different shaft segments on the shafting and the waveguide transfer relationship between the connected shaft segments are established based on the waveguide theory, and the waveguide coefficient is estimated combined with the vibration response observation data. The estimation of longitudinal oil film stiffness of thrust bearing and admittance of thrust bearing housing and hull connection structure is avoided, and the method of propeller longitudinal excitation based on waveguide transfer relation is verified by a specific example. The results show that the method is feasible and effective, and provides a new method for accurately grasping the actual propeller excitation.


Introduction
Due to the limitations of the working environment and structural conditions, the propeller excitation of the marine propulsion shafting is difficult for direct measurement [1]. The spectrum of the unsteady force of the propeller is composed of the spectrum of the low-frequency discrete force with the blade-passing-frequency (BPF) and double-BPF generated by the uneven force on the blades, and the spectrum of the low-frequency broadband force generated by the interaction between the stern turbulence and the propeller, which can be obtained mainly through experimental tests, semi-empirical formulas [2], and numerical calculations [3]. In the test pool, the scaled model can be measured to obtain the spectrum of the unsteady force of the propeller, but the correlation between the test results and the real ship model is still unclear. Besides, the lateral force measurement is more complicated, and no effective measurement means have been established yet. The numerical calculation of the spectrum of the low-frequency discrete force mainly includes the lifting surface method, panel method, and CFD method, etc. [4]. In the lifting surface method, the blade is assumed as a thin wing to calculate the pressure distribution of the blade, but the influence of the hub is not taken into consideration. In the panel method, the calculation of the propeller excitation force can reflect the actual operating conditions of the propeller to the greatest extent, the influence of the hub is taken into account, more accurate calculation results can be obtained, and there is no need for the thin wing hypothesis [5,6]. Based on the software Fluent for fluid mechanics computation, Ma Chao applied the RNG turbulence model to analyze the uneven turbulent flow near the propeller and the non-uniform pressure pulsation near the propeller under the influence of the turbulent flow. The author also studied the variation of the bearing force and moment of the propeller [7]. It is difficult to calculate the low-frequency broadband force spectrum of the propeller. Some software for fluid computation can be used to theoretically simulate and analyze the propeller excitation in the flow field, or building a physical model of the shafting in the laboratory to conduct experimental research on the propeller excitation characteristics. Still, they are different from the real shaft system. Therefore, grasping the propeller excitation under actual working conditions is very important for maneuvering the vibration and noise of the propeller and propulsion shafting, and controlling the vibration of the hull structure caused by the shafting vibration.
Based on the waveguide theory, the waveguide model of different shaft sections of the shafting and the waveguide transfer relationship between the connected shaft sections are established in response to the above difficulties. Assuming there is a shaft section suitable for direct vibration response measurement. The waveguide coefficient can be estimated based on the vibration response observation data of the shaft section and the pre-established shaft section waveguide model. The longitudinal excitation of the propeller can be derived from the waveguide model of the shaft section and each connecting shaft section of the propeller and the waveguide transfer relationship between these connecting shaft sections. Since the waveguide model and the vibration response observation data of the shaft section are used to estimate the waveguide coefficients, the boundary conditions of the shaft section are not involved, so the estimation of the longitudinal oil film stiffness of the thrust bearing and the admittance of the joint structure between the thrust bearing housing and the hull are avoided. Figure 1 is a general multi-shaft (rod) section longitudinal vibration model based on the general structure of the shafting, wherein the propulsion shafting is assumed to consist of N uniform shaft sections, and x0, x1, x2, …, xN represent coordinates of the ednpoint of the shaft (rod) sections. Note that the length of the i-th (i=1, 2, …, N) shaft (rod) section is li, and the cross-sectional area is Si. Each rod section has a complex elastic modulus of Ei * =Ei(1+ji) and a density of ρi. The mass of the propeller is mP. The coupling stiffness between the shafting and the hull (thrust bearing axial stiffness) is kA. The admittance in the longitudinal vibration direction of the joint position between the shafting and the hull is YA. The longitudinal excitation of the propeller is Fa. Let u(x, t) represent the longitudinal vibration displacement of the particle on the shaft, and it satisfies the following wave function and boundary conditions:

Theoretical Model Construction
Let ui(x, t) represent the vibration displacement of each shaft section, then where Ui represents the Fourier transform of ui, and the section stiffness is given by ( ) The waveguide method [8] is adopted to establish the shaft model (as shown in Figure 1), and each shaft section is described by wave solution Ui(x, ): where i is the complex wavenumber, and i 2 =ρiω 2 /Ei * . At the far left end (set x0=0), the boundary conditions are given by: Therefore, if C1 + and C1 − are measurable or identifiable, the longitudinal excitation Fa of the propeller can be calculated by the above equation. Now suppose that u1(x, t) or U1(x, ) of the first shaft section at the left side can be measured arbitrarily, then arrange several measuring points on the shaft at the location of x=x1 (1) , x2 (1) ,x3 (1) , …, xn (1) , we get: where, j (1) presents the measurement error. These n measuring points are presented as a matrix expression ( ) Applying the principle of least square estimation, the residual expression is Using the above equation to estimate C1 + and C1 − requires n2. When the shaft section is simplified to a uniform cross-section structure and the model parameters are accurate, Eqs. (5), (7), and (10) can be used to calculate the longitudinal excitation of the propeller. Because theoretical derivation is rigorous, there is no need for verification with specific examples. However, since the real shafting system usually has many structures such as key grooves, local bumps, and external coating materials, various degrees of simplification are needed in the simplifying process of the shaft section into an ideal homogeneous, uniform cross-section waveguide model, which will introduce theoretical calculation errors. In addition, the measurement of vibration response will also generate errors. To reduce the impact of measurement errors, the least square estimation algorithm of Eq. (10) is adopted.

The Problem Wherein Measuring Points Cannot Be Arranged on The Shaft Section Close to The Propeller
The method of using Eqs. (5), (7), and (10) to estimate the longitudinal excitation of the propeller from the measured longitudinal vibration response does not require the parameters of the other connecting shaft sections after the first shaft section in Figure 1, and the longitudinal stiffness of the thrust bearing of the right side of the shafting and housing admittance parameters, which is the convenience of the above method. However, in practical applications, it may be challenging to arrange sensors on the stern section of the shafting close to the propeller to measure the longitudinal vibration response. Therefore, it is necessary to further consider picking up the longitudinal vibration response on the shaft section close to the thrust bearing to verify the feasibility of the indirect estimation of the longitudinal excitation of the propeller.
It is still assumed that each shaft section is a constant-section structure, and its longitudinal  (3). At the junction of any i-th section and its left (i-1)-th section (i2), there are continuity conditions for displacement response and internal force, namely Therefore, if the Ci + and Ci − of the i-th shaft section are measurable, then the Ci−1 + and Ci−1 − of the (i-1)-th shaft section are calculated from the above formula. Using this recurrence relationship, C1 + and C1 − can be derived. Finally substitute C1 + and C1 − into Eq. (5) to estimate the longitudinal excitation Fa.
For the measurement of Ci + and Ci − , the least square estimation method of Eq. (10) is still used, namely, where U (i) =[U1 (i) , U2 (i) ,U3 (i) , …, Un (i) ] T is a column vector composed of the longitudinal vibration response at n2 measuring points arranged on the i-th shaft section. Refer to Eq. (7), R (i) can be written as It should be pointed out that the application of the above recursive algorithm requires as much as possible to ensure the accuracy of the waveguide model of the measurable shaft section and that on its left side up to the propeller. In actual operation, approximating parameters in the waveguide model will cause calculation errors in each recursive calculation, and the errors will be accumulated during multiple recursions. Therefore, the arrangement of measuring points on the shaft far away from the propeller will increase the uncertainty in theoretical calculations and reduce the estimation accuracy of the propeller excitation. Although the recursive algorithm is theoretically feasible, the measuring points should be arranged as close as possible to the longitudinal excitation shaft section.

Measuring Waveguide Coefficients of Any Shaft Section by Combining Waveguide Method and Modal Method
For the forward problem of dynamics, the complete solution of the shaft waveguide model requires clear boundary conditions and external excitation [9]. However, the waveguide coefficients of the shaft section can be estimated through the measured response, so the boundary conditions of the waveguide model are not necessary, which allows unknown kA at the right end of the shafting ( Figure  1) and admittance YA in the longitudinal vibration direction at the connection between the thrust bearing housing and the hull. The basic principle of combining the waveguide method and the modal method to measure the waveguide coefficients of any shaft section is to identify the influence factors of the main vibration modes through the arrangement of as few measuring points as possible, then apply the modal superposition theory to fit the analytical expression of Ui(x, ) curve, and finally calculate Ui(x, )/x on this basis. The detailed steps are as follows: Firstly, the shaft section to be measured is a certain part of the entire shafting. Separating the shaft section in the entire shafting, that is, releasing the constraints of other shaft sections or structures in the shafting to the shaft section to be measured and replacing them with binding forces FL and FR. The shaft section to be measured is regarded as a free rod at both ends, and FL and FR are regarded as the external excitation. Even the longitudinal vibration mode of this shaft section is not accessible for theoretically calculation; its numerical solution can be easily obtained through accurate finite element modeling, that is, when the shaft section to be measured is regarded as a free rod at both ends, its longitudinal vibration mode is easy to be determined with guaranteed accuracy.
Let j (i) (x) represent the mode function of each order of the longitudinal vibration of the measured shaft section. By using the modal superposition theory, we get where, Mj (i) and j (i) represent the modal quality and modal frequency of the shaft section, respectively; i is the damping loss factor; j (i) represents the impact factor (weight factor) of the j-th mode when the modal superposition is used to synthesize Ui.
where 1, 2, …, P represent the modal truncation error generated when Ux . In fact, it is the core of the problem to estimate the low-frequency components in the longitudinal excitation of the propeller, and the natural frequencies of the shaft sections are generally relatively high. For example, for a 1-m steel shaft with a uniform cross-section, a free boundary at both ends, except for the modal frequency of its rigid body f0=0, the lowest non-zero natural frequency is f12594.4 Hz. Therefore, even the impact factors of the first two modes (presumably designing the first three low-frequency modes, including the rigid body mode) in the low-frequency band are included, it can be assured the modal truncation error is within a smaller order of magnitude. In the real test analysis, the main component of the error vector  will be the measurement error.
In order to eliminate the influence of the error vector  as much as possible, the least square estimation principle is applied to Eq. (18) to estimate the modal influence factor, and we have The modal influence factor can be estimated with Eq. (19). Generally, multiple observation points should be arranged in a scattered manner, and the number is larger than that of main influencing modes, i.e., P>Q. The more observation points there are, the higher the confidence in the estimation of . For estimating the low-frequency longitudinal excitation of the propeller, as mentioned above, it is usually sufficient to consider the first three low-frequency modes including the rigid body mode, that is, the number of observation points P>3, which is relatively easy to achieve in the real test. Besides, care should be taken to avoid placing the observation point on (or near) the node of the main influencing mode, because this will lead to the impact factor of this mode cannot be truthfully reflected in the observation data.
After the impact factor of the primary influencing mode is estimated, the response function and the slope function of the shaft section to be measured can be expressed as where ( ) i j  represents the estimated value of the impact factor (the j-th element of ξ ) of the j-th primary influencing mode.
Substitute the above equation into Eq. (20), the we obtain the estimated value of the waveguide coefficients as where g1 and g2 are differential operators, with Finally, it must be pointed out that by the convergence theorem of Fourier series [10], Eq. (21) cannot be used to estimate the waveguide coefficient near the end of the rod; To do this, you should first estimate waveguide coefficient at the middle of the rod using Eq. (21), and then calculate the waveguide coefficient at the end with Eq. (14).

Examples
The following calculation examples are used to theoretically verify the method of estimating the waveguide model through modal superposition method.
Assume a round steel shaft (rod) with a uniform cross-section and both free ends with the length L=1 m, radius r=0.03 m, complex elastic modulus E * =21010 9 (1+0.01j) Pa, density ρ=7800 kg/m 3 . Establish a one-dimensional coordinate system along the rod axis, and set x=0 at the left end and x=L at the right end. Excitation Extend the above equation into the following even function with a period of 2L And then expand ( ) , Ux  to a Fourier series, for 0xL, we have Note that substituting Eqs. (27) and (28) The first term of the above equation corresponds to the impact factor of the longitudinal vibration rigid body mode (natural frequency 0=0, mode type 0=1) of a uniform rod with free ends.
According to Eq. (29), the (complex) natural frequency of the rest non-rigid body modes, mode type, and mode impact factor respectively are Therefore, for a uniform free rod excited from both ends, the wave solution (Eq. (24)) is an exact analytical solution, and the modal superposition solution (Eq. (28)) is the Fourier series expansion of the wave solution. If the excitation frequencies of FL and FR are low (generally lower than the first-order non-zero natural frequency of the rod), because the high-order mode has a small impact factor, the superposition synthesis with finite low-order modes is used to approximate U(x, ). Figures 2 and 3 show the comparison between the U(x, ) by superposition synthesis with the first Q-order modes (including rigid body modes) and its precise wave solution. It is assumed that FL=cost (N) and FR=−sint (N). The excitation frequency of f=/(2)=100 Hz, 1000 Hz is selected.  Figure 2. Comparison of U(x, ) by superposition synthesis with the first low-order modes (including rigid body mode) and their precise wave solutions (excitation frequency f=100 Hz).
Based on Figures 2 and 3, in the low-frequency band, the consistency of U(x, ) by superposing several low-frequency modes and the precise wave solution is mainly determined by the mode order involved in modal superposition, but slightly affected by the change of excitation frequency. Obviously, the more modal orders involved in the modal superposition, the closer the result of modal superposition is to the exact wave solution. Besides, notice that at the left and right ends of the rod, there are significant deviations between the finite modal superposition result and the exact wave solution (the point of x=0 in (a) and (c), and x=L in (b) and (d) in Figures 2, 3 and 4). It can be attributed to the usually discontinuous derivatives at x=0 and x=L, obtained by extending U(x, ) to get ( ) , Ux  through Eq. (25). In fact, from the Fourier series solution of Eq. (29), there will be U(x, )/x0 at x=0 and x=L, but the wave solution of Eq. (24) does not lead to the same conclusion, which can be explained by referring to the convergence theorem of the Fourier series.
In summary, the accuracy of U(x, ) synthesized by superposing finite-order modes is mainly determined by the convergence of the Fourier series, which has a slower convergence speed at both