Estimation of longitudinal excitation of propeller using a novel hybrid method

1. Introduction

The propeller excitation force is an essential exciting source for submarine underwater navigation. When the propeller runs in a temporally varying and spatially uneven flow field near the stern of the submarine, the hull can be excited to produce noise radiation by the unsteady dynamic force caused by the non-uniform force on each blade, which has an important impact on the underwater acoustic stealth performance of the submarine [1]. How to simulate the real stern flow conditions and the study of the scale effect and similarity between the experimental model and the real model are essential to accurately predicting the direct noise of the propeller.

Due to the limitations of the working environment and structural conditions, the propeller excitation of the marine propulsion shafting is difficult for direct measurement [2]. 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 [3], and numerical calculations [4]. 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. [5]. 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 [6, 7]. 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 [8]. 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.

According to the dynamics theory, the excitation can also be estimated by measuring the vibration response of the system when the real excitation cannot be measured, which belongs to the second type of inverse problem of the dynamics [9]. Since it is difficult to directly measure the fluid excitation of the propeller in working conditions in seawater, it becomes a reasonable choice to estimate the propeller excitation by measuring the vibration response of the shafting. This estimation must be carried out on the premise that the vibration model of the shafting has been established. Therefore, a theoretical model for estimating the excitation of the propeller from the shaft vibration response should be first established based on the general structure of the propulsion shafting. However, there exist some practical difficulties in dealing with the above problem solely by the modal theory, for examples:

1. It is not easy to arrange a large number of vibration sensors on the entire shafting, especially the stern section close to the propeller. Therefore, it is almost impossible to fully grasp the actual vibration of the shafting through the vibration response measurement.

2. Although it is technically feasible to establish a finite element model (FEM) based on the structural size and material properties of the shafting, the determination of the boundary conditions of the shafting is difficult, especially for the shaft-hull connections with thrust bearings. There is currently a lack of reliable theoretical calculations or measurement methods for the parameters, including the longitudinal oil film stiffness $k_A$ of the thrust bearing and the admittance $Y_A$ of the joint structure between the thrust bearing housing and the hull.

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. At the same time, after the waveguide model of the entire shafting is established, the method for estimating the waveguide coefficient of a shaft section using the vibration response can also be used to calculate the longitudinal excitation from the shafting to the hull through the thrust bearing, which is the same as the principle of estimating the longitudinal excitation of the propeller. Based on this, the admittance function test problem of thrust bearing housing can also be discussed.

2. Basic principle

2.1. Theoretical model construction

The unsteady force generated by the propeller in the non-uniform wake field at the stern of the hull, the radial simple harmonic component force generated by the main engine device and the longitudinal component force generated by the torsional vibration of the propulsion shafting all make the propulsion shafting produce periodic longitudinal tension and compression deformation, which is the common longitudinal vibration of the propulsion shafting of large ships and submarines. The analysis of the longitudinal vibration mechanism of propulsion shafting focuses on solving two problems: one is to calculate the longitudinal vibration characteristics of propulsion shafting and its influencing factors, and to judge whether there is longitudinal resonance of propulsion shafting within the normal speed range of the main engine; The second is to analyze the longitudinal dynamic response of the typical excitation of the propulsion shafting, and judge whether the response amplitude is within the allowable range, so as to avoid the adverse effects of the secondary excitation force generated by the longitudinal vibration of the shafting on other structures of the hull.

Fig. 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 $x_0$, $x_1$, $x_2$, …, $x_N$ represent coordinates of the endpoint of the shaft (rod) sections. Note that the length of the $i$-th ($i=$1, 2, …, $N$) shaft (rod) section is $l_i$, and the cross-sectional area is $S_i$. Each rod section has a complex elastic modulus of $E_i^{*}=E_i\left(1+j{\eta }_i\right)$ and a density of ${\rho }_i$. The mass of the propeller is $m_P$. The coupling stiffness between the shafting and the hull (thrust bearing axial stiffness) is $k_A$. The admittance in the longitudinal vibration direction of the joint position between the shafting and the hull is $Y_A$. The longitudinal excitation of the propeller is $F_a$.

Fig. 1Simplified model of the longitudinal vibration of the propulsion shafting

Let $u\left(x,t\right)$ represent the longitudinal vibration displacement of the particle on the shaft, and it satisfies the following wave function and boundary conditions:

1

$\rho S\left(x\right)\frac{{\partial }^{2}u}{\partial t^2}-\frac{\partial }{\partial x}\left[E^{*}\left(x\right)\cdot S\left(x\right)\frac{\partial u}{\partial x}\right]=F_a\delta \left(x\right).$

By measuring the longitudinal vibration response of the shafting, the longitudinal excitation of the propeller can be estimated mainly based on the following considerations:

(1) Lack of direct measurement means for the longitudinal excitation $F_a$ of the propeller;

(2) For real systems, the connection stiffness $k_A$ is generally not easy to be determined;

(3) The relevant geometric parameters and material characteristic parameters of each shaft section and the propeller of the propulsion shafting can be accurately obtained;

(4) The longitudinal vibration response of the shaft system can be directly measured.

Let $u_i\left(x,t\right)$ represent the vibration displacement of each shaft section, then:

2

${\rho }_i{S}_i{\ddot{u}}_i-E_i^{\mathrm{*}}{S}_i{u}_{i,xx}=0\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{o}\mathrm{r}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}{\omega }^2{\rho }_i{S}_i{U}_i\left(x,\omega \right)+E_i^{\mathrm{*}}{S}_i{U}_{i,xx}\left(x,\omega \right)=0,$
$\left({x_i}_{-1}\le x\le x_i,i=1,2,\cdots ,N\right),$

where $U_i$ represents the Fourier transform of $u_i\text{,}$ and the section stiffness is given by $\overline{\left(ES\right)}=\sum _{k=1}^{\kappa }\overline{E_k^{*}}\overline{S_k}$.

The waveguide method [10] is adopted to establish the shaft model (as shown in Fig. 1), and each shaft section is described by wave solution $U_i\left(x,\omega \right)$:

3

$U_i\left(x,\omega \right)=C_i^{+}\left(\omega \right)\mathrm{e}\mathrm{x}\mathrm{p}\left(-j{\lambda }_{i}x\right)+C_i^{-}\left(\omega \right)\mathrm{e}\mathrm{x}\mathrm{p}\left(j{\lambda }_{i}x\right),\left({x_i}_{-1}\le x\le x_i,\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}i=1,2,\cdots ,N\right),$

where ${\lambda }_i$ is the complex wavenumber, and ${l_i}^2={\rho }_i{\omega }^2/{E_i}^{\mathrm{*}}$. At the far left end (set $x_0=$0), the boundary conditions are given by:

4

$m_P{\ddot{u}}_1\left(x_0,t\right)=F_a+E_1^{*}{S}_1{u}_{1,x}\left(x_0,t\right)\mathrm{}\mathrm{o}\mathrm{r}-{\omega }^2{m}_P{U}_1\left(x_0,\omega \right)=F_{\mathrm{a}}+E_1^{\mathrm{*}}{S}_1{U}_{1,x}\left(x_0,\omega \right).$

Substitute Eq. (3) into Eq. (4), then we have:

5

$F_a=-{\omega }^2{m}_P{U}_1\left(x_0,\omega \right)-E_1^{*}{S}_1{U}_{1,x}\left(x_0,\omega \right)$
$=-{\omega }^2{m}_P\left[C_1^{+}\left(\omega \right)+C_1^{-}\left(\omega \right)\right]-j{\lambda }_1{E}_1^{*}{S}_1\left[-C_1^{+}\left(\omega \right)+C_1^{-}\left(\omega \right)\right]$
$=C_1^{+}\left(\omega \right)\left(-{\omega }^2{m}_P+j{\lambda }_1{E}_1^{*}{S}_1\right)-C_1^{-}\left(\omega \right)\left({\omega }^2{m}_P+j{\lambda }_1{E}_1^{*}{S}_1\right).$

Therefore, if $C_1^{+}$and $C_1^{-}$ are measurable or identifiable, the longitudinal excitation $F_a$ of the propeller can be calculated by the above equation. Now suppose that $u_1\left(x,t\right)$ or $U_1\left(x,\omega \right)$ 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=x_1^{\left(1\right)}$, $x_2^{\left(1\right)}$, $x_3^{\left(1\right)}$, …, $x_n^{\left(1\right)}$, we get:

6

$U_1\left(x_j^{\left(1\right)},\omega \right)=C_1^{+}\left(\omega \right)\mathrm{e}\mathrm{x}\mathrm{p}\left(-j{\lambda }_1{x}_j^{\left(1\right)}\right)+C_1^{-}\left(\omega \right)\mathrm{e}\mathrm{x}\mathrm{p}\left(j{\lambda }_1{x}_j^{\left(1\right)}\right)+{\epsilon }_j^{\left(1\right)},j=\mathrm{1,2},\dots ,n,$

where, ${\epsilon }_j^{\left(1\right)}$ presents the measurement error. These $n$ measuring points are presented as a matrix expression:

7

${\mathbf{U}}^{\left(1\right)}=\left[\begin{array}{l}{U}_1^{\left(1\right)}\\ U_2^{\left(1\right)}\\ ⋮\\ U_n^{\left(1\right)}\end{array}\right]=\left[\begin{array}{cc}\mathrm{e}\mathrm{x}\mathrm{p}\left(-j{\lambda }_1{x}_1^{\left(1\right)}\right)& \mathrm{e}\mathrm{x}\mathrm{p}\left(j{\lambda }_1{x}_1^{\left(1\right)}\right)\\ \mathrm{e}\mathrm{x}\mathrm{p}\left(-j{\lambda }_1{x}_2^{\left(1\right)}\right)& \mathrm{e}\mathrm{x}\mathrm{p}\left(j{\lambda }_1{x}_2^{\left(1\right)}\right)\\ ⋮& ⋮\\ \mathrm{e}\mathrm{x}\mathrm{p}\left(-j{\lambda }_1{x}_n^{\left(1\right)}\right)& \mathrm{e}\mathrm{x}\mathrm{p}\left(j{\lambda }_1{x}_n^{\left(1\right)}\right)\end{array}\right]\left[\begin{array}{l}{C}_1^{+}\left(\omega \right)\\ C_1^{-}\left(\omega \right)\end{array}\right]+\left[\begin{array}{c}{\epsilon }_1^{\left(1\right)}\\ {\epsilon }_2^{\left(1\right)}\\ ⋮\\ {\epsilon }_n^{\left(1\right)}\end{array}\right]={\mathbf{R}}^{\left(1\right)}{\mathbf{C}}_1+{\mathbf{\epsilon }}^{\left(1\right)}.$

Applying the principle of least square estimation, the residual expression is:

8

${\left[{\mathbf{\epsilon }}^{\left(1\right)}\right]}^T{\mathbf{\epsilon }}^{\left(1\right)}={\left[{\mathbf{U}}^{\left(1\right)}-{\mathbf{R}}^{\left(1\right)}{\mathbf{C}}_1\right]}^T\left[{\mathbf{U}}^{\left(1\right)}-{\mathbf{R}}^{\left(1\right)}{\mathbf{C}}_1\right].$

Let:

9

$\nabla \left\{{\left[{\mathbf{\epsilon }}^{\left(1\right)}\right]}^T{\mathbf{\epsilon }}^{\left(1\right)}\right\}=\nabla \left\{{\left[{\mathbf{U}}^{\left(1\right)}-{\mathbf{R}}^{\left(1\right)}{\mathbf{C}}_1\right]}^T\left[{\mathbf{U}}^{\left(1\right)}-{\mathbf{R}}^{\left(1\right)}{\mathbf{C}}_1\right]\right\}$
$=2\nabla \left\{{\left[{\mathbf{U}}^{\left(1\right)}-{\mathbf{R}}^{\left(1\right)}{\mathbf{C}}_1\right]}^T\right\}\left[{\mathbf{U}}^{\left(1\right)}-{\mathbf{R}}^{\left(1\right)}{\mathbf{C}}_1\right]=-2{\left[{\mathbf{R}}^{\left(1\right)}\right]}^T\left[{\mathbf{U}}^{\left(1\right)}-{\mathbf{R}}^{\left(1\right)}{\mathbf{C}}_1\right]=0,$

we have:

10

${\mathbf{C}}_1={\left\{{\left[{\mathbf{R}}^{\left(1\right)}\right]}^T{\mathbf{R}}^{\left(1\right)}\right\}}^{-1}{\left[{\mathbf{R}}^{\left(1\right)}\right]}^T{\mathbf{U}}^{\left(1\right)}.$

Using the above equation to estimate $C_1^{+}$and $C_1^{-}$ requires $n\ge$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.

3. Measuring waveguide coefficients of any shaft section by combining waveguide method and modal method

The general longitudinal vibration displacement of the shafting is described by Eq. (3). Its calculation accuracy depends entirely on the calculation accuracy of the section stiffness $\overline{\left(ES\right)}$ and wavenumber $\lambda$, that is, the accuracy of the physical model of the shafting. As for determining the waveguide coefficients of $C^{+}$and $C^{-}$, it can boil down to the initial value problem of the differential equation of Eq. (4). As long as the accuracy of the initial value and the physical model of the shafting is sufficient, the obtained waveguide coefficients can be guaranteed to have enough calculation accuracy. Lastly, Eq. (14) describes the conversion relationship of the waveguide coefficients when each shaft section of the shafting has a uniform cross-section. The calculation accuracy also completely depends on the physical modeling accuracy of the shafting. Therefore, the longitudinal excitation of the propeller can be estimated using the above-mentioned waveguide transfer model if the waveguide coefficients of some section of the shafting are given.

For the forward problem of dynamics, the complete solution of the shaft waveguide model requires clear boundary conditions and external excitation [11]. However, to calculate the propeller excitation problem from the measured shaft vibration response, it must be assumed that the vibration response of a certain section of the shafting can be measured, and 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 $k_A$ at the right end of the shafting (Fig. 1) and admittance $Y_A$ 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 $U_i\left(x,\omega \right)$ curve, and finally calculate $\partial U_i\left(x,\omega \right)/\partial 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 $F_L$ and $F_R$. The shaft section to be measured is regarded as a free rod at both ends, and $F_L$ and $F_R$ 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 ${\mathrm{\Phi }}_j^{\left(i\right)}\left(x\right)$ 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, $M_j^{\left(i\right)}$ and ${\omega }_j^{\left(i\right)}$ represent the modal quality and modal frequency of the shaft section, respectively; ${\eta }_i$ is the damping loss factor; ${\xi }_j^{\left(i\right)}$ represents the impact factor (weight factor) of the $j$-th mode when the modal superposition is used to synthesize $U_i$.

Suppose that $P$ observation points of $x_1^{\left(i\right)}$, $x_2^{\left(i\right)}$, …, $x_P^{\left(i\right)}$are taken from the shaft section. The measured values of longitudinal vibration response are ${\tilde{U}}_i\left(x_1^{\left(i\right)},\omega \right)$, ${\tilde{U}}_i\left(x_2^{\left(i\right)},\omega \right)$, …, ${\tilde{U}}_i\left(x_{}^{\left(i\right)},\omega \right)$, and according to Eq. (17) the measured value of each observation point can be expressed as the following matrix equation:

where ${\epsilon }_1$, ${\epsilon }_2$, …, ${\epsilon }_P$ represent the modal truncation error generated when ${\tilde{U}}_i\left(x_1^{\left(i\right)},\omega \right)\text{,}$${\tilde{U}}_i\left(x_2^{\left(i\right)},\omega \right)$, …, ${\tilde{U}}_i\left(x_{}^{\left(i\right)},\omega \right)$ are synthesized with the $Q$-order low-frequency modes, and the measurement error in the actual observation value ${\tilde{U}}_i\left(x_1^{\left(i\right)},\omega \right)$, ${\tilde{U}}_i\left(x_2^{\left(i\right)},\omega \right)$, …, ${\tilde{U}}_i\left(x_{}^{\left(i\right)},\omega \right)$. 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 $f_0=$0, the lowest non-zero natural frequency is $f_1\approx$ 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 $\epsilon$ will be the measurement error.

In order to eliminate the influence of the error vector $\epsilon$ 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 $\mathrm{\epsilon }$. 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 ${\tilde{\xi }}_j^{\left(i\right)}$ represents the estimated value of the impact factor (the $j$-th element of $\tilde{\xi }$) 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 $g_1$ and $g_2$ are differential operators, with:

Finally, it must be pointed out that by the convergence theorem of Fourier series [12], 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).

4. 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⁹(1+0.01j) Pa, density $\rho =$7800 kg/m³. 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.

4.1. The relationship between wave solution and modal superposition solution

Excitation $F_L\left(\omega \right)$ and $F_R\left(\omega \right)$ are applied onto the left and right end of the rod. The longitudinal vibration displacement $U\left(x,\omega \right)$ of the rod satisfies the waveguide model of Eq. (2). Because rod cross-section is uniform, the waveguide coefficients $C^{+}\left(\omega \right)$ and $C^{-}\left(\omega \right)$ are independent of $x$. According to the following boundary conditions:

22

$E^{*}S{\left.\left(\frac{\partial U}{\partial x}\right)\right|}_{x=0}=F_L,E^{*}S{\left.\left(\frac{\partial U}{\partial x}\right)\right|}_{x=L}=F_R,$

the waveguide coefficient can be obtained as:

23

$C^{+}=\frac{F_R-F_L\mathrm{e}\mathrm{x}\mathrm{p}\left(j\lambda L\right)}{j\lambda E^{*}S\left[-\mathrm{e}\mathrm{x}\mathrm{p}\left(-j\lambda L\right)+\mathrm{e}\mathrm{x}\mathrm{p}\left(j\lambda L\right)\right]}，C^{-}=\frac{F_R-F_L\mathrm{e}\mathrm{x}\mathrm{p}\left(-j\lambda L\right)}{j\lambda E^{*}S\left[-\mathrm{e}\mathrm{x}\mathrm{p}\left(-j\lambda L\right)+\mathrm{e}\mathrm{x}\mathrm{p}\left(j\lambda L\right)\right]}.$

Substitute the above equation into Eq. (3), we get:

24

$U\left(x,\omega \right)=\frac{\left[F_R-F_L\mathrm{e}\mathrm{x}\mathrm{p}\left(j\lambda L\right)\right]\mathrm{e}\mathrm{x}\mathrm{p}\left(-j\lambda x\right)+\left[F_R-F_L\mathrm{e}\mathrm{x}\mathrm{p}\left(-j\lambda L\right)\right]\mathrm{e}\mathrm{x}\mathrm{p}\left(j\lambda x\right)}{j\lambda E^{*}S\left[-\mathrm{e}\mathrm{x}\mathrm{p}\left(-j\lambda L\right)+\mathrm{e}\mathrm{x}\mathrm{p}\left(j\lambda L\right)\right]}.$

Extend the above equation into the following even function with a period of 2$L$:

25

$\overline{U}\left(x,\omega \right)=\left\{\begin{array}{l}\begin{array}{ll}U\left(x,\omega \right),& \left(0\le x\le L\right),\end{array}\\ \begin{array}{ll}-U\left(x,\omega \right),& \left(-L\le x\le 0\right),\end{array}\end{array}\right.$
$\overline{U}\left(x+2L,\omega \right)=\overline{U}\left(x,\omega \right).$

And then expand $\overline{U}\left(x,\omega \right)$ to a Fourier series, for $0\le x\le L$, we have:

26

$U\left(x,\omega \right)=\frac{a_0}{2}+\sum _{j=1}^{\infty }{a}_j\mathrm{c}\mathrm{o}\mathrm{s}\frac{j\pi x}{L},$

where:

27

$a_0=\frac{2}{L}{\int }_0^{L}U\left(x,\omega \right)dx$
$={\int }_0^L\frac{\left[F_R-F_L\mathrm{e}\mathrm{x}\mathrm{p}\left(j\lambda L\right)\right]\mathrm{e}\mathrm{x}\mathrm{p}\left(-j\lambda x\right)+\left[F_R-F_L\mathrm{e}\mathrm{x}\mathrm{p}\left(-j\lambda L\right)\right]\mathrm{e}\mathrm{x}\mathrm{p}\left(j\lambda x\right)}{j0.5\lambda E^{*}SL\left[-\mathrm{e}\mathrm{x}\mathrm{p}\left(-j\lambda L\right)+\mathrm{e}\mathrm{x}\mathrm{p}\left(j\lambda L\right)\right]}dx$
$=\frac{-\frac{\left[F_R-F_{L}exp\left(j\lambda L\right)\right]}{j\lambda }{\left.\mathrm{e}\mathrm{x}\mathrm{p}\left(-j\lambda x\right)\right|}_0^L+\frac{\left[F_R-F_L\mathrm{e}\mathrm{x}\mathrm{p}\left(-j\lambda L\right)\right]}{j\lambda }{\left.\mathrm{e}\mathrm{x}\mathrm{p}\left(j\lambda x\right)\right|}_0^L}{j0.5\lambda E^{*}SL\left[-\mathrm{e}\mathrm{x}\mathrm{p}\left(-j\lambda L\right)+\mathrm{e}\mathrm{x}\mathrm{p}\left(j\lambda L\right)\right]}$
$=2\frac{-F_L+F_R}{{\left(j\lambda \right)}^2{E}^{*}SL}=2\frac{F_L-F_R}{\rho SL{\omega }^2},$

28

$a_j=\frac{2}{L}{\int }_0^{L}U\left(x,\omega \right)\mathrm{c}\mathrm{o}\mathrm{s}\frac{j\pi x}{L}dx=\frac{\left[F_R-F_L\mathrm{e}\mathrm{x}\mathrm{p}\left(j\lambda L\right)\right]{\int }_0^L\mathrm{e}\mathrm{x}\mathrm{p}\left(-j\lambda x\right)\mathrm{c}\mathrm{o}\mathrm{s}\frac{j\pi x}{L}dx}{j0.5\lambda E^{*}SL\left[-\mathrm{e}\mathrm{x}\mathrm{p}\left(-j\lambda L\right)+\mathrm{e}\mathrm{x}\mathrm{p}\left(j\lambda L\right)\right]}$
$+\frac{\left[F_R-F_L\mathrm{e}\mathrm{x}\mathrm{p}\left(-j\lambda L\right)\right]{\int }_0^L\mathrm{e}\mathrm{x}\mathrm{p}\left(j\lambda x\right)\mathrm{c}\mathrm{o}\mathrm{s}\frac{j\pi x}{L}dx}{j0.5\lambda E^{*}SL\left[-\mathrm{e}\mathrm{x}\mathrm{p}\left(-j\lambda L\right)+\mathrm{e}\mathrm{x}\mathrm{p}\left(j\lambda L\right)\right]}=\frac{2}{L}\frac{F_R\mathrm{c}\mathrm{o}\mathrm{s}\left(j\pi \right)-F_L}{E^{*}S\left[{\left(j\lambda \right)}^2+{\left(\frac{j\pi }{L}\right)}^2\right]}$
$=\frac{F_R\mathrm{c}\mathrm{o}\mathrm{s}\left(j\pi \right)-F_{\mathrm{L}}}{0.5\rho SL\left[-{\omega }^2+{\left(\frac{j\pi }{L}\right)}^2\frac{E^{*}}{\rho }\right]}.$

Note that substituting Eqs. (27) and (28) into Eq. (26), we get the modal superposition equation, namely:

29

$U\left(x,\omega \right)=-\frac{F_{\mathrm{R}}-F_{\mathrm{L}}}{\rho SL{\omega }^2}+\sum _{j=1}^{\infty }\frac{F_{\mathrm{R}}\mathrm{c}\mathrm{o}\mathrm{s}\left(j\pi \right)-F_{\mathrm{L}}}{0.5\rho SL\left[-{\omega }^2+{\left(\frac{j\pi }{L}\right)}^2\frac{E^{*}}{\rho }\right]}\mathrm{c}\mathrm{o}\mathrm{s}\frac{j\pi x}{L}.$

The first term of the above equation corresponds to the impact factor of the longitudinal vibration rigid body mode (natural frequency ${\omega }_0=$0, mode type ${\mathrm{\Phi }}_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:

30

${\omega }_j=\frac{j\pi }{L}\sqrt{\frac{E^{*}}{\rho }},{\mathrm{\Phi }}_j\left(x\right)=\mathrm{c}\mathrm{o}\mathrm{s}\frac{j\pi x}{L},{\xi }_j\left(\omega \right)=\frac{F_{\mathrm{R}}\mathrm{c}\mathrm{o}\mathrm{s}\left(j\pi \right)-F_{\mathrm{L}}}{0.5\rho SL\left[-{\omega }^2+{\left(\frac{j\pi }{L}\right)}^2\frac{E^{*}}{\rho }\right]}.$

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 $F_L$ and $F_R$ 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\left(x,\omega \right)$. Figs. 2-4 show the comparison between the $U\left(x,\omega \right)$ by superposition synthesis with the first $Q$-order modes (including rigid body modes) and its precise wave solution. It is assumed that $F_L=\mathrm{c}\mathrm{o}\mathrm{s}\omega t$ (N) and $F_R=-\mathrm{s}\mathrm{i}\mathrm{n}\omega t$ (N). The excitation frequency of $f=\omega /2\pi =$100 Hz, 400 Hz, 1000 Hz is selected.

Fig. 2Comparison of Ux,ω by superposition synthesis with the first low-order modes (including rigid body mode) and their precise wave solutions (excitation frequency f= 100 Hz)

a) Real part of $U\left(x,\omega \right)$

b) Imaginary part of $U\left(x,\omega \right)$

c) Real part of $U\left(x,\omega \right)$

d) Imaginary part of $U\left(x,\omega \right)$

Fig. 3Comparison of Ux,ω by modal superposition (MS) synthesis of the first low-order modes (including rigid body mode) and their precise wave solutions (excitation frequency f= 400 Hz)

a) Imaginary part of $U\left(x,\omega \right)$

b) Imaginary part of $U\left(x,\omega \right)$

c) Real part of $U\left(x,\omega \right)$

d) Imaginary part of $U\left(x,\omega \right)$

Based on the Figs. 2-4, in the low-frequency band, the consistency of $U\left(x,\omega \right)$ 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 Figs. 2-4). It can be attributed to the usually discontinuous derivatives at $x=$0 and $x=L$, obtained by extending $U\left(x,\omega \right)$ to get $\overline{U}\left(x,\omega \right)$ through Eq. (25). In fact, from the Fourier series solution of Eq. (29), there will be $\partial U\left(x,\omega \right)/\partial x\equiv$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\left(x,\omega \right)$ synthesized by superposing finite-order modes is mainly determined by the convergence of the Fourier series, which has a slower convergence speed at both ends of the rod.

Fig. 4Comparison of Ux,ω by superposition synthesis with the first low-order modes (including rigid body mode) and their precise wave solutions (excitation frequency f= 1000 Hz)

a) Real part of $U\left(x,\omega \right)$

b) Imaginary part of $U\left(x,\omega \right)$

c) Real part of $U\left(x,\omega \right)$

d) Imaginary part of $U\left(x,\omega \right)$

4.2. Estimating waveguide coefficients based on the vibration response data of finite observation points and modal function

Note that the $U\left(x,\omega \right)$ obtained by the modal superposition method (Eq. (29)) cannot be used to calculate the end excitation $F_L$ and $F_R$ of the rod, because it does not satisfy Eq. (22). Therefore, in order to use the vibration response data of the limited observation points to invert the end excitation of the rod, the waveguide coefficients $C^{+}$and $C^{-}$ must be first estimated, and then substituted into the wave solution of Eq. (3); finally, $F_L$ and $F_R$ can be calculated using Eq. (22). In addition, for the uniform cross-section rod in this example, Eqs. (7)-(10) can be used to estimate the waveguide coefficients. This method does not need to be studied through simulation examples. The following is mainly based on Eqs. (18)-(21) to estimate the waveguide coefficient by combining the waveguide method and the modal method.

Fig. 5The influence of the number of mode superposition on the estimation accuracy of waveguide coefficients (Excitation frequency f= 100 Hz)

a) Real part of $U\left(x,\omega \right)$, at $P=$ 4, $Q=$ 3

b) Imaginary part of $U\left(x,\omega \right)$, at $P=$ 4, $Q=$ 3

c) Real part of $U\left(x,\omega \right)$, at $P=$ 5, $Q=$ 4

d) Imaginary part of $U\left(x,\omega \right)$, at $P=$ 5, $Q=$ 4

e) Real part of $U\left(x,\omega \right)$, at $P=$ 7, $Q=$ 6

f) Imaginary part of $U\left(x,\omega \right)$, at $P=$ 7, $Q=$ 6

Firstly, as in Eq. (18), suppose that the first $Q$-order modes of the rod are used to synthesize $U\left(x,\omega \right)$, and arrange $P$ ($P\ge Q$) observation points on the rod. Since the modal superposition series solution of $U\left(x,\omega \right)$ converges slowly at the rod end, it is necessary to avoid arranging observation points close to the end of the rod; Temporarily select a section of 2/3 rod length from the middle of the rod for arranging the measurement points, and arrange the P observation points evenly.

Secondly, substitute the coordinate value of each observation point into Eq. (24), and calculate the longitudinal vibration response of each observation point as the theoretical simulation of the observation value to form the observation vector $\tilde{\mathbf{U}}$ in Eq. (18). Substitute the coordinate value of each observation point into the mode function of Eq. (30), and then take the first $Q$-order mode function value to form the observation point mode type matrix $\mathrm{\Phi }$ of Eq. (18), and then estimate the mode impact factor vector $\tilde{\mathbf{\xi }}$ by Eq. (19). Further substitute $\tilde{\mathbf{\xi }}$ and the coordinate values of each observation point into Eq. (21) to estimate the waveguide coefficients ${\tilde{C}}^{+}$ and ${\tilde{C}}^{-}$. Note that because rod cross-section is a uniform cross-section, the ${\tilde{C}}^{+}$ and ${\tilde{C}}^{-}$ should be constant in theory. But due to the calculation errors, the calculated waveguide coefficients corresponding to different observation points are different in practice, so the mean value of the estimated waveguide coefficients will be used as the final estimated value of ${\tilde{C}}^{+}$ and ${\tilde{C}}^{-}$.

Finally, in the following Figs. 5-7, the exact wave solution of $U\left(x,\omega \right)$, $Q$-order modal superposition estimation and wave solution estimation based on the $Q$-order modal superposition estimation are compared to illuminate the feasibility of the proposed method for estimating the waveguide coefficients and the factors affecting the accuracy of the estimation results. The exact wave solution of $U\left(x,\omega \right)$ is calculated by Eq. (24). The $Q$-order modal superposition estimation of $U\left(x,\omega \right)$ is obtained by substituting the modal impact factor vector $\tilde{\mathbf{\xi }}$ into the first equation of Eq. (20). The wave solution estimation based on the $Q$-order modal superposition estimation is calculated by substituting the estimated value of the waveguide coefficient into Eq. (3).

Fig. 6The influence of the change of excitation frequency (in the low-frequency domain) on the accuracy of waveguide coefficient estimation

a) Real part of $U\left(x,\omega \right)$, at $f=$ 400 Hz

b) Real part of $U\left(x,\omega \right)$, at $f=$ 1000 Hz

Fig. 5 shows a set of figures, indicating that the more modal orders included in the modal superposition, the higher the accuracy of the estimation of the waveguide mode; but as the modal order increases, its influence in the low-frequency domain decreases rapidly, thus for low-frequency problems, there is no need for too many modal orders. Since the required number $P$ of observation points is no less than the order $Q$ of the mode participating in the superposition, and $P=Q+1$ is used in Figs. 5(a)-(e). Fig. 6 shows that when the excitation frequency changes in the low-frequency band, the estimation accuracy of the waveguide model is marginally affected. Fig. 7 shows that when the number of modal orders is small, the estimation accuracy of the waveguide coefficient can also be increased by increasing the number of measuring points (Note that no matter the number of measuring points, you should avoid placing measuring points near the end of the rod). However, when the modal order is sufficiently large, the number $P$ of measurement points ($P\ge Q$) has no significant effect on the estimation accuracy of the waveguide coefficient.

Fig. 7The influence of the number of observation points on the estimation accuracy of waveguide coefficients (Excitation frequency f= 1000 Hz)

a) Real part of $U\left(x,\omega \right)$, at $P=$ 5, $Q=$ 3

b) Real part of $U\left(x,\omega \right)$, at $P=$ 7, $Q=$ 3

c) Real part of $U\left(x,\omega \right)$, at $P=$ 6, $Q=$ 4

d) Imaginary part of $U\left(x,\omega \right)$, at $P=$ 8, $Q=$ 4

5. Conclusions

The waveguide method is used to estimate the excitation of the propeller. Because there could be many intermediate transmission shaft sections between the measuring shaft section and the propeller, it is necessary to accurately describe the waveguide model of each shaft section and the transfer relationship between them. Based on the multi-section shaft model of Fig. 1, the obtained conclusions are summarized as follows:

1) If the shaft section connected with the propeller can be regarded as a uniform shaft section, and its longitudinal vibration response can be measured, the waveguide coefficients $C_1^{+}$ and $C_1^{-}$ can be simply estimated with Eqs. (7)-(10), and then the longitudinal excitation of the propeller can be calculated by using Eq. (5); Because the measuring point is close to the excitation position, this estimation displays high reliability.

2) In practice, when it is challenging to arrange sensors close to the stern section near the propeller, it is necessary to establish the waveguide model and the waveguide coefficient transfer relationship between the measuring shaft section and each shaft section connected to the propeller. And the propeller longitudinal excitation can be calculated by using Eqs. (5) and (14)-(16).

3) In multiple recursive calculations, the approximate processing of the parameters of the waveguide model can lead to the accumulation of calculation errors. Therefore, arranging measurement points far away from the shaft section near the propeller can reduce the estimation accuracy of the propeller excitation, and the measurement points should be arranged as close as possible to the shaft section near the propeller.

4) This work is limited to the theoretical demonstration of a hybrid method for estimating the longitudinal propeller excitation from the measured longitudinal vibration response. Because the theoretical derivation is rigorous, numerical simulation to verify the aforementioned method is no longer necessary, and its practical implementation is also feasible. However, designing a specific experimental system to conduct certain experimental research to clarify the specific influencing factors for errors in the actual excitation estimation is still valuable.