System Identification and Dynamic Analysis of the Propulsion Shaft Systems Using Response Surface Optimization Technique

1. Introduction

The propulsion shaft system is the main component of a marine power system consisting of a motor, thrust, intermediate and stern bearings, couplers, and propeller. Its primary use is to provide thrust and propulsion forces to marine vehicles to maneuver [1]. The propulsion shaft system experiences dynamic vibration excitations during operations, leading to longitudinal, lateral, and torsional vibrations in the shaft [2]. Axial excitation is caused by the non-uniform flow of the stern near the propeller in time and space, resulting from the asymmetry of the hull, which leads to a longitudinal vibration of the shaft [2]. Similarly, lateral vibration in the system is irrefutable and is induced by various reasons, such as the rotation of the propeller in the non-uniform wake that causes lateral excitation in the shaft. In addition, the lateral vibration in the shaft is also caused by the unbalanced force transmitted to the hull via bearings and misalignment between components such as bearings, shafts, and couplers [3]. Likewise, torsional vibration in the marine propulsion shaft is induced by the cyclic torque of the power system, which breaks the shaft, resulting in the rupture of couplers and bolts [2].

In the past, many researchers have attempted different approaches to model the dynamics of marine propulsion shaft systems. The two most common methods to model the dynamics are Transfer Matrix Method (TMM) and Finite Element Method (FEM). TMM models the dynamics of the structure by dividing the structure into smaller segments, and the relationship between the smaller segments is represented using transfer matrices. These transfer matrices are then solved subsequently to obtain system dynamics. TMM can be used in two ways for the analysis of dynamic systems: transfer matrix relating to the state vector [4] and transfer matrix relating to the constant coefficients of differential equation solutions [5]. The advantage of the latter, as compared to the transfer matrix method relating to state vector when applied to shaft structure, is the possibility of reducing the number of multiplied matrices when adjacent shaft segments have the same material properties and diameters. Chahr-Eddine and Yassine [5] used TMM related to the constant coefficients of differential equation solutions to study the force axial and torsional vibrations of a shaft line. In their study, they investigated the normal and tangential stress tensor components resulting from axial-torsional deformations and vibrations in the propeller and intermediate shafts under the influence of propeller-induced static and variable hydrodynamic excitations. Their result [5] shows that the strength of the shaft line depends on the value of the static tangential stresses. More about the TMM approach to model the dynamic of the marine shaft can be found in references [6,7]. The Finite Element Method (FEM) involves discretizing the complex structure in the space dimensions to grids called finite elements through the construction of a mesh. Each element represents a portion of the system with a finite number of degrees of freedom at nodes. The response of the entire system is approximated by assembling these elements without overlapping, like constructing complex structures from simpler components. Li et al. [6] developed the dynamic model of axial vibration of a coupled propeller shaft system using FEM. Similarly, Zhang et al. [8] used FEM and TMM to study the longitudinal vibration in the marine propulsion shaft system structure, where the author discussed the variation in first natural frequency due to bearings stiffness, length of the shaft, and number of bearings. Huang et al. [9] developed a FEM of marine propulsion shaft systems that incorporates a coupled constraint on the propeller elements. The authors studied the torsional and transverse vibration of the system under idling and loading conditions at different rotational speeds. Their research provides insight into the basic principles of marine propulsion shafting coupled dynamics. Their study also supports the prediction of coupled vibration, which can improve the safety and reliability of ship sailing performance. Yucel and Arpaci [10] analyzed the ship hull structural vibration by creating a three-dimensional FEM, which includes the ship hull, deckhouse, and machinery propulsion system. The FEM model was used for local and global vibration analyses under free-free (dry) and in-water (wet) conditions. The wet analysis utilized acoustic elements. To account for overall ship hull structure vibration, a combination of several damping components was considered for total damping in their study. More on the dynamics studies in the propeller shaft system can be found in the literature [11,14,17].

Upon conducting a thorough review of the existing literature, it is found that while FEA models are extensively employed to study the dynamics of propulsion shaft systems, there is a significant drawback in utilizing the FEA to construct dynamic models. In particular, the model response can deviate significantly when FE results are not correlated with experimental results. A marine propeller shaft system essentially consists of several components, including bearings, couplers, etc., whose stiffness values are unknown. These unknown stiffness parameters will alter the system’s dynamics entirely if the FE model is not correlated with experimental results using a suitable optimization technique. This research is focused on identifying the unknown dynamics parameters (stiffness) of the system via a suitable optimization algorithm that minimizes the error between the measured and FEA responses. The response surface optimization (RSO) technique is proposed in this paper. RSO is a mathematical optimization technique that has been utilized in the past for optimizing process parameters in casting and welding, as well as reliability and fatigue optimization studies of electronic packages [12,13,18,19]. However, the application of RSO in a complex system such as a marine shafting line is still in its infancy. The RSO algorithm implemented in this paper is integrated into the ANSYS Workbench environment and is readily accessible, optimizing simulation data to match existing experimental data using parametric optimization. RSO offers several advantages over other optimization techniques, such as direct optimization [20], including computational time efficiency, parameter space exploration, optimization flexibility, and model interpretability.

The rest of the paper is organized as follows. Section 2 presents the structural model updating using the RSO algorithm and the validation of the RSO algorithm in a 3-degree of freedom (DOF) system. The experimental configuration and approach of generating lateral and longitudinal response in a propeller shaft system is presented in Section 3. Section 4 presents the FE model development of the propeller shaft system and implementation of the RSO algorithm using the experimental data. Finally, Section 5 concludes the paper with the findings of the research work.

2. Structural Model Updating Using Response Surface Optimization

The basic idea of the response surface-based optimization technique utilized in this paper is to update the initial FE model with an estimated model. The key stages are as follows: (i) selecting updating parameters, potentially utilizing sensitivity analysis; (ii) sampling updating parameters via the design of experiment (DOE) technique and computing the response using the FE model; (iii) creating a response surface through regression analysis between the updating parameters and the associated response, accompanied by regression error analysis; (iv) building objective functions based on the simulated and measured response of the structure; and (v) repeatedly iterating and optimizing objective functions within the established response surface model.

The selected updating parameters should be able to clarify the ambiguity of the model. If the number of updating parameters exceeds the number of structural responses available, an ill-conditioned optimization problem may appear. To remove the ill-conditioned, sensitivity analysis of the model parameters can be done and only crucial parameters can be selected to update the model. Furthermore, the sampling of updating parameters affects the accuracy and computation efficiency of the response surface model. For this, a commonly used approach called Design of Experiment (DOE) is used. Various methodologies are available to calculate DOE, such as CCD (Central Composite Design), Box-Behnken Design, Sparse Grid Initialization, etc. In this research, the Latin Hypercube Sampling (LHS) method with Full Quadratic Model Samples is chosen as the DOE method, which is solved to obtain the value of output parameters. The response surface generated by genetic aggregation can be expressed as an ensemble utilizing the weighted average of different metamodels [13].

$${\widehat{y}}_{ens}\left(x\right)=\sum _{i=1}^N{w}_i.\widehat{y_i}\left(x\right)$$

(i)

Where, ${\widehat{y}}_{ens}$ is the prediction of the ensemble, $\widehat{y_i^{}}$ is the prediction of the i^th response surface,$N$ is the number of metamodels used, and $w_i$is the weight factor of the i^th response surface.

The weight factors satisfy the following criteria:

$$\sum _{i=1}^N{w}_i = 1 and w_i\ge 0,1\le i\le N$$

(ii)

The best weight factor is estimated by minimizing the root mean square error (RMSE) between the actual and the predicted values. A cross-validation is utilized by taking the predicted residual error sum of squares (PRESS), which is calculated for a number of candidate models for the same data set, with the lowest value of PRESS indicating the best model. The RMSE and PRESS can be calculated by using equation Equations (iii) and (iv).

$$RMSE{(\widehat{y}}_{ens})=\sqrt{\frac{1}{N}\sum _{j=1}^N(y\left(x_j\right)-{\widehat{y}}_{ens}{\left(x_j\right))}^2}$$

(iii)

$$PRES{S}_{RMSE}{(\widehat{y}}_{ens})=\sqrt{\frac{1}{N}\sum _{j=1}^N(y\left(x_j\right)-{\widehat{y}}_{ens,-j}{\left(x_j\right))}^2}$$

(iv)

$${\widehat{y}}_{ens,-j}\left(x\right)=\sum _{i=1}^N{w}_i.{\widehat{y}}_{i,-j}\left(x\right)$$

(v)

Where, x_j is the j^th design point, $y\left(xj\right)$is the output parameter value at$x_j,{\widehat{y}}_{i,-j}$is the prediction of the i^th response surface without the j^th design point and $N$ is the number of design points for the design of experiment.

2.1. Numerical Example to Demonstrate the Response Surface Optimization (RSO) Technique

2.1.1. Response Data Generation and Implementation of RSO

To demonstrate the working principle of response surface optimization, we comprehensively analyze a 3 Degree of Freedom (DOF) discrete mass system model with lumped parameters. We use two approaches: forward and reverse. The forward approach involves creating the model’s initial response, while the reversed approach is used to determine the model’s parameters based on response data. To simulate the model’s initial response, we create a forced response model in the ANSYS harmonic analysis module with predetermined values for the spring stiffness. In the reversed approach, we consider the spring stiffness as an unknown parameter and calculate it by analyzing the model’s response, including the maximum amplitude response and resonant frequencies. To optimize the unknown parameters, we use the Latin hypercube sampling technique and generate response surfaces using the genetic aggregation technique. Finally, we compare the optimized parameters to the initial assumed stiffness based on the least RMSE and PRESS, confirming the effectiveness of the RSO technique. Figure 1 illustrates the spring mass model used in the analysis.

The equation of the Motion for a system shown in Figure 1 can be written as,

$$\left[m\right]\left\{\ddot{x}\right\}+\left[k\right]\left\{x\right\}=\left\{F\right\}$$

(vi)

Where $\left[m\right]$and $\left[k\right]$are the mass and stiffness matrix, $\left\{F\right\}$ is the external excitation force vector. If the damping matrix $\left[C\right]$ is introduced in the model defined by Equation (vi), then Equation (vi) can be written as,

$$\left[m\right]\left\{\ddot{x}\right\}+\left[k\right]\left\{x\right\}+\left[C\right]\dot{\left\{x\right\}}=\left\{F\right\}$$

(vii)

Equation (vii) can be solved by using the FE based harmonic analysis tools utilizing the frequency domain approach [15,16] or by using the numerical technique. In this paper ANSYS modal and harmonic analysis module are used to generate the response data. The parameters used in the simulation are as follows.

$$m_1=1kg,m_2=2kg,m_3=1kg,k_1=5.7\times {10}^5\frac{N}{m},k_2=4.3\times {10}^5\frac{N}{m},$$

$k_3=7.5\times {10}^5\frac{N}{m},k_4=8.5\times {10}^5\frac{N}{m}$ and modal damping ratios ${\zeta }_1={\zeta }_2={\zeta }_3=0.05$ are assumed for the simulation.

The resonant frequencies obtained from the Modal analysis are listed in Table 1.

The harmonic analysis of a system shown in Figure 1 is carried out by exciting$m_3$ with a force of $50N$ harmonic at an excitation frequency range of (50-250) Hz (covering all three resonant frequencies), and the steady state amplitude at each frequency is recorded. The response of Mass 2 at the excitation frequency range is shown in Figure 2 as a Frequency Response Function (FRF). Having known one of the FRFs (Figure 2) and the resonant frequency shown in Table 1, the reverse approach for parameter estimation is carried out by assuming the four stiffness of the spring mass system to be unknown.

2.1.2. Reverse Approach for Parameter Identification

In the reverse analysis approach, the spring-mass system’s stiffness parameters are assumed to be unknown. The simulated response data from the forward approach shown in Figure 2 and the resonant frequencies presented in Table 1 are used as the output to obtain the unknown stiffness parameters by creating the response surfaces. The outputs utilized are the first mode, second mode, third mode resonant frequencies and the maximum acceleration responses at the corresponding resonant frequencies. The goal is to obtain the unknown stiffness parameters of the spring through the design of experiments such that the difference between the measured output and the reconstructed output is minimal. For this, the design of experiment (DOE) is created using 4 spring stiffness parameters, which are the dependent variable, and the upper limit and lower limit of all the spring’s stiffness are assigned as $1\times {10}^{5}to9.5\times {10}^{5}N/m$. To establish the correlation between the stiffness of each of the four springs and the desired output, 25 design points are generated to create the response surfaces. After the optimization, three candidate models are obtained; they are then cross verified using the PRESS. The candidate model with the lowest value of PRESS is taken as the best model. The associated stiffness parameters of the best model and the comparison with the original spring stiffness is tabulated in Table 2. As shown in Table 2, the error between the original stiffness and the optimized stiffness is less than 2%, indicating the robustness of the implemented optimization technique. Figure 3, Figure 4 and Figure 5) show the comparison of the harmonic response of the model from the forward approach and the reversed approach. It should be noted that only mass 2 response data and 3 resonant frequencies are taken for optimization. The results will be identical if mass 1 or mass 3 response data are taken for analysis. Similarly, Figure 6 and Figure 7 show the response plot between the dependent variable and the objective function. While the spring stiffness variation is taken from $1\times {10}^5$ to $9\times {10}^{5}N/m$ during DOE, the wide range of variation can be taken even starting from 100 to$1\times {10}^{9}N/m$; the only difference will be the computational cost; one needs to create a huge number of design points, resulting in a longer time for convergence.

3. Experimental Analysis of Propulsion Shaft System

An experimental setup was designed to demonstrate the system identification concept using the response surface optimization algorithm. The proposed experimental model was fabricated in-house by downscaling the suitable available model in the existing literature [17]. Figure 7 shows the test model of the propulsion system used for longitudinal excitation testing. The details of the experimental setup consist of the propeller, intermediate shafts, thrust bearing, radial contact ball bearings, different types of couplings, and the variable frequency drive motor. The variable frequency drive motor was utilized during the experiment to support the shaft. The structure’s right side has an intermediate shaft connected to the motor by a jaw coupler. A thrust bearing and the first intermediate bearing follow it. Then, a rigid coupler connects the first intermediate shaft to the second one. Finally, a propeller is attached to the end of the second intermediate shaft using set screws. The material properties and the detailed dimensions of the individual parts are listed in Table 3. To carry out forced vibration tests, we utilized a modal exciter from PCB Piezotronics. This exciter has a peak sine force of 58 N (13 lbf) and a peak-peak stroke of 18 mm (0.7”). It also features an adjustable collet stinger attachment and a through-hole armature design, as shown in Figure 7. We attached an ICP impedance head to the exciter for the measurement of the input acceleration/force excitation signal. Experiments were conducted in two arrangements: (i) lateral or transverse excitation, and (ii) longitudinal or axial excitation.

3.1. Transverse (Lateral) Vibration Excitation

Lateral vibration in the propulsion shaft system occurs due to misalignment of the shaft and manufacturing defects in the gearing system. In order to mimic the lateral vibration in the laboratory environment, an experimental setup was designed by using an exciter that directly excites the shaft in the transverse direction. Figure 8 shows our experimental design for the transverse vibration setup, where the intermediate shaft was excited using the modal shaker, resulting in the lateral vibration of the entire assembly. Two PCB (352C22) accelerometers were used to record the transverse vibration of the shaft: one behind the propeller and the other exactly at the intermediate shaft face, as shown in Figure 8.

During the experiment, we utilized harmonic excitation within the frequency range of 50-1000 Hz, incorporating a 1 Hz sweep sine input and an acceleration excitation amplitude of 0.1 g, implementing an open loop controller. The steady state response of two accelerometers was recorded at each excitation frequency using the Crystal Instruments EDM software. To identify the vibration mode of the assembly, we plotted the response of the shaft system within the given excitation frequency range using the steady state response of the shaft at each excitation frequency. Figure 9 displays the lateral response of the shaft located near the propeller and away from the propeller. As anticipated, the response of the shaft close to the propeller was greater than that of the shaft located away from the propeller, as there was no support system in place near the propeller.

3.2. Longitudinal Vibration Excitation

For the longitudinal experimental setup, the propeller end was excited with very low excitation with an amplitude of 0.1 g to imitate the real-world excitation scenario provided by the unsymmetrical stern flow near the propeller due to the current in the water. Two accelerometers were attached near the propeller hub face, and another one was attached to the face of the rigid coupling, 450 mm away from the propeller. The purpose of this setup was to analyze the structural characteristics in the vicinity of the propeller and to observe any changes in vibration response along the shaft. Figure shows the vibration response of the shaft near and away from the propeller at an excitation acceleration of 0.1 g. Since the shaft is short as compared to the practical application, it is expected that the longitudinal mode of the vibration of shaft is high as compared to the practical scenario, which is evident from Figure 11.

4. Finite Element Modeling of Propulsion Shaft System

The FEM of the propulsion shaft system is carried out in ANSYS Release 23 Workbench Environment. The goal is to define the suitable boundary conditions and connections in the assembly while creating the finite element model such that the structure’s dynamic behavior is as similar as possible between the actual tested model response and the numerically simulated response. The detailed 3D model of the propulsion shaft system is modeled using SolidWorks and is imported to the ANSYS Workbench environment. The imported simplified model for FEM analysis is shown in Figure 11. As illustrated in Figure 11, the model consists of two ball bearing housing that connects with the shaft using four springs each, and the thrust bearing is suppressed where four additional springs are added as a body-to-ground connection in a shaft. The boundary conditions imposed are three fixed supports: one on the end of the shaft, which connects to the motor; the other two fixed supports are used on the faces of the bearing housing where the screws are used to fix it to the support as shown in the experimental model in Figure 7 and Figure 8 respectively.

Selection of Bearing Stiffness Parameters and Updating the Parameters by using Response Surface Optimization

The FEM model updating process starts by choosing a suitable set of numerical model parameters (unknown bearings stiffness) that are efficiently changed during the process. These parameters represent unknown properties in the experiment. The fewer unknown parameters, the more efficient the FEM. If the number of unknown parameters exceeds the output parameters, the optimization may produce irrelevant and ill-conditioned results. Some unknown parameters in a system have a more significant impact on a structure’s properties than others. Parameter correlation looks at how sensitive each input parameter is to the output parameter, which is essential in optimizing a model. Parameters with low correlation can be given less priority or ignored to reduce the computational time needed during optimization. The model that we consider consists of 12 springs. The first four springs represent the modeled thrust ball bearing, and the connection type used is the body to ground. The other eight spring connections are connected from body to body from the shaft to the bearing housing (Two Ball Bearings). This is done to represent the contact between the inner ring of the bearing and the shaft as well as to create a parameter for the parametric study to find the unknown longitudinal (ANSYS uses the term Longitudinal stiffness for all the radial stiffness) stiffness, which is then used as Design of Experiments (DOE) to create response surface and then lead to optimization. These 12 longitudinal stiffness are used as input parameters and provided with lower and upper bounds, where the lower bound is set to ${10}^5$ N/m and the upper bound is set to ${10}^9$ N/m. A total of 3 output parameters are chosen: the shaft’s first mode resonant frequency, the second mode resonant frequency, and its maximum amplitude response obtained from the transverse vibration test. Since only 3 output parameters are utilized to update the structural model, and there are 12 unknown stiffnesss, the problem is ill-conditioned. It is well known that not all the bearing stiffness in four directions [9,17] for each bearing will contribute to the transverse vibration in the excitation frequency range (50 -1000 Hz) of interest; for this, a design of experiment is carried out in the ANSYS optimization module. The fundamental procedure in the design of experiment consists of determining projected values of the response features at various design points (DP) in the parameter space by simulating each DP. This response surface will be a placeholder for the FE simulation model (Figure 12). The number of DP created depends upon the response surface method to be established. The DOE of test sample points is used to fit the response surface, which is the point in response surface technique. A good selection of sample points can reduce the calculation cost of DOE and improve the accuracy of the response surface. The commonly used point selection method at this stage is to fill the design space with the most effective and minimum sample points as much as possible. The positions of the test sample points meet certain symmetry and uniformity requirements. We utilize the Latin hypercube sampling method with full quadratic sample methods to create the design points through which the response surface is generated. Once the response surface is generated, local sensitivity analysis is conducted to see how each input parameter (bearing stiffness) affects the output response. This information helps to eliminate parameters that do not correlate with the output. Figure 13 shows the sensitivity analysis of the 12 stiffness parameters, of which only three stiffness parameters are sensitive by more than 20%. So, an ill-conditioned problem is synthesized into a well-posed problem by optimizing those three spring parameters.

After synthesizing the problem into a well-posed one, the optimization process involves defining objectives and constraints for the output parameters. The target value is selected from the experimental measurements obtained, as shown in Figure 9. To carry out the optimization, the MOGA method (Multi-Objective Genetic Algorithm) is used, which is a variant of the popular NSGA-II (Non-dominated Sorted Genetic Algorithm-II) based on controlled elitism concepts. This method supports multiple objectives and constraints and aims to find the global optimum. Figure 14 displays the response of the optimized model and the experimental model. The error between the two models is less than 10% in terms of the resonant frequency and the maximum amplitude of vibration, which is considered an acceptable range for a complex problem with several unknown variables. It’s worth noting that only three stiffness parameters of the bearing contributed to the transverse vibration problems in the chosen excitation frequency range. The other stiffness parameters may contribute to the longitudinal excitation or other higher modes that are not excited.

5. Conclusions

This paper presents the detailed dynamic analysis and system identification technique of a marine propulsion shaft system subjected to dynamic excitation. The detailed dynamic analysis in a propulsion shaft system is conducted by designing a suitable scaled experimental setup that can induce significant vibration when subjected to dynamic excitations. The system identification technique is carried out using the response surface optimization technique. First, a numerical example is presented to demonstrate the concept of response surface optimization using the 3 DoF lumped parameter model. Following the numerical example, the same algorithm is utilized to update the FE model in a propulsion shaft system using the experimental data measured through harmonic excitation. The results show that the ill-conditioned problem can be resolved to a well-posed problem by performing the sensitivity analysis, and the unknown parameters of the structure can be estimated by using the design of experiment techniques in the FE-based environment. The efficiency of this approach lies in its ability to resolve complex problems and generate valuable insights that can drive innovation and progress.