Coupling characteristic analysis of ship shafting design parameters and research on multidisciplinary design optimization

3.1. Whirling vibration sub-disciplinary model

The rotor dynamics equation of propeller propulsion shafting is:

where, [$M$] is the mass matrix, considering the water-attached mass of propeller, [$C$] is the damping matrix, [$K$] is the stiffness matrix, {$F$} is the force acting on the shafting, and [$G$] is the gyro effect matrix. When {$F$} is a zero matrix, the eigenvalue of the homogeneous equation system is the natural frequency of the shafting and the solution is the free vibration response. The forced vibration response of shafting can be calculated by applying the approximate propeller force {$F$}.

Based on the simplification and analysis above, the corresponding finite element model of shafting is established on Ansys platform, as shown in Fig. 5. Beam188 is used to establish the Axial section and Combin14 is used to establish the support spring system. Bearing 214 is adopted to represent the longitudinal spring element at thrust bearing. The elastic modulus of shafting material $E$, 210 GPa, Poisson’s ratio $\mu$, and density $\rho$ is 210 GPa. 0.32 and 7850 kg/m³ respectively. The buoyancy coefficient of the immersion section $\kappa$ is –1000 kg/m³.

Fig. 5Shafting finite element model

The external input of whirling vibration calculation is propeller bearing force. According to the theory of propeller lifting line, the propeller is divided into $Nm$ segments along the radial direction. When the blade angle is $\theta \text{'}$, the thrust $F_x$, torque $M_x$ and tangential force $F_{{\theta }^{\text{'}}}$ of the blade can be expressed as:

where, $r_c$ is radial position of segment $c$, $C$ is the number of segments in radial direction, $Z$ is the number of propeller blades. The amplitude of propeller bearing force at each rotation speed was obtained by quasi-steady calculation and applied to the shafting. According to instruction of shafting components and shafting structural parameters in Table 1, the shafting geometry model was established based on CATIA platform. The approximate calculation of propeller thrust and radial force was carried out in MATLAB. The calculation process of cyclotron vibration is shown in Fig. 6. Each bearing model is considered as a spring with vertical and lateral stiffness. Fixed restraint is imposed at the thrust bearing end, welding constraints are imposed on the thrust bearing end, and spring restraint is applied at the bearing. By applying the calculation force at the propeller end, the whirling vibration was calculated by the harmonic vibration module. Considering the vibration characteristics of the long shaft system, all analyzes are implemented within 0-200 Hz.

Fig. 6Simulation flow of whirling vibration

The whirling vibration response of each bearing monitoring point under initial state is obtained from the numerical results of shafting vibration model, as shown in Fig. 7.

Fig. 7Bearing amplitude-frequency characteristics of whirling vibration at initial state

It can be seen from Fig. 7 that, in the range of 0-60 Hz, the amplitude of whirling vibration of the rear stern bearing is much higher than that of the front stern bearing and thrust bearing, and the amplitude of $Y$ direction is higher than that of $Z$ direction. Therefore, great priority should be given to the improvement of whirling vibration characteristics of the stern bearing in the optimization process.

3.2. Alignment sub-discipline model

Each shaft segment is discretized into a combination of elastic beam elements in the alignment calculation, for the shafting is considered as a multi-support continuous beam with equal cross-section, ignoring the radius change of the shaft section. Each beam element consists of two nodes to connect adjacent elements.

Fig. 8Stress analysis diagram of beam element

According to the relevant theory of material mechanics, for any element a, as shown in Fig. 8, $J$ and $K$ is the node in the left and right end surfaces of the element respectively, and then the displacement and force of the element can be expressed as:

where, ${\omega }_j$, ${\omega }_k$, ${\theta }_j$, ${\theta }_k$ represent the linear displacement and angular displacement of the node in the direction of y-axis deflection, respectively; $T_j$, $M_j$, $T_k$, $M_k$ represent the shear force and bending moment of the node $j$ and $k$, respectively; and ${\left[K\right]}^e$is the element stiffness matrix:

where, $\beta$ is the influence coefficient of shear deformation, $A$ and $G$ are shear modulus of elasticity; then, the relationship $f$ between shafting force and deformation can be obtained, where $K$ is the stiffness matrix of shafting, which can be obtained by combining $L$. After obtaining the stiffness matrix, the actual geometric parameters of shafting and the support stiffness can be calculated. The alignment analysis of the shafting in the initial state (linear alignment) is carried out. The calculation results are shown in Fig. 9.

From Fig. 9, it can be seen that the shafting stern shaft deflection is the largest, with the maximum value of 0.479 mm. The maximum stress occurs at the cross-section change part of the stern shaft, with the maximum value of 15.93 MPa. The rear stern bearing load (62.19 kN) is higher than the front stern bearing load (18.25 kN) and the thrust bearing load (25.12 kN). It shows that the load and vibration characteristics of the rear stern bearing are poor, which should be given priority in the optimization process.

Fig. 9Calibration results of shafting alignment in initial state

a) Rotation angle

b) Deflection

c) Equivalent stress

3.3. Sub-disciplinary model of bearing oil film stiffness

Bearing support stiffness is an important parameter in shafting alignment calculation and vibration calculation. It can be expressed as series values of static stiffness and oil film dynamic stiffness [25]:

where, $K_{sum}$ is total equivalent stiffness at bearing, $K_{const}$ is stiffness of basic structure, $K_{oil}$ is stiffness of oil film. Ship shafting usually operates at low speed and heavy load. Based on hydrodynamic lubrication theory, the Reynolds equation can be used to approximate the dynamic characteristics of oil film. The working principle of radial lubrication bearing in ship shafting is shown in Fig. 10.

The oil film distributes in the gap between the rotor and the bearing bush when the rotor runs in the bearing. The oil film fluid can be considered as the laminar flow applicable to the Reynolds equation, so the liquid film equation between the shaft diameter and the shaft bush is:

Assuming that the oil film pressure at the bearing boundary is 0 and the shaft section has noinclination inside the bearing, the radial force $F_{bd}$ and tangential force $F_{bt}$ of the bearing under stable operation can be expressed as follows:

where, $R_b$ is the bearing radius, $\omega$ is the shaft speed, $\mu$ is the lubricant oil viscosity coefficient, $L_b$ is the bearing length, $\stackrel{-}{r}$ is the average radial clearance inside the bearing, $\epsilon$ is the bearing eccentricity.

The resultant force $F_{oil}$ could be calculated from the radial force and tangential force of oil film acts on the bearing through the bearing bush, which can be expressed as:

The load of the shafting can be obtained by the calibration calculation model established above. The magnitude and direction of the load are the same as that of $F_{oil}$ numerically. Further derivation can obtain the expressions of bearing eccentricity ɛ and force attitude angle $\stackrel{-}{\theta }$:

where, $S_f$ is modified Sommerfeld number. From Eq. (13), the remaining parameters can be approximately solved according to any two parameters of load, rotation speed and eccentricity. In order to simplify the calculation and construct a stiffness model of oil film support suitable for multidisciplinary optimization design, the bearing stiffness is simplified to a linear spring in both horizontal and vertical directions, and the stiffness matrix of oil film can be approximately expressed as follows:

The lateral stiffness $K_{zz}$ and the vertical stiffness $K_{yy}$ can be expressed as:

where $h$ is the approximate value of oil film thickness.

3.4. Decoupling approximate response analysis of design parameters

The coordinated operation of the shafting model requires high calculation costs. In order to accelerate the optimization speed, uniform design is adopted for the design variables. Taking the calculation results of each sub-model as sample points, the radial basis function neural network is used to generate the target response surface. Neural network has been widely used in many fields for its advantage of cost saving [26]. In present study, 5000 initial sample points are generated in the design variable space for training and verification. The load-design variable response surface of rear stern bearing is obtained as shown in Fig. 11.

Fig. 11Force response surface of rear stern bearing

a) Force response surface of $d$ and $y_1$

b) Force response surface of $d$ and $y_2$

c) Force response surface of $d$ and $y_3$

d) Force response surface of $y_1$ and $y_2$

e) Force response surface of $y_1$ and $y_3$

f) Force response surface of $y_2$ and $y_3$

From Fig. 11, it is shown that the coupling response surface of the shaft inner diameter $D$ and the bearing vertical position $Y$ is nearly linear, while the coupling response surface between the vertical positions of the shafts is non-linear. The load of the tail bearing decreases with the increase of the radius. The higher the vertical position of the bearing, the greater the load of the stern bearing, indicating that the thrust bearing has the greatest influence coefficient on the whole shafting. When the thrust bearing is too high, the intermediate bearing will be empty.

Fig. 12Stiffness response surface of rear stern bearing oil film

a)$Y$-direction stiffness response surface of rear stern bearing

b)$Z$-direction stiffness response surface of rear stern bearing

According to the approximate equation of oil film stiffness characteristics, it can be seen that the physical and chemical characteristics of lubricating oil itself are relatively stable under stable operation, and the stiffness characteristics of oil film are mainly affected by shaft rotation speed and bearing load. Secondary approximate response surface of bearing load-speed-oil film stiffness is further obtained by RBF neural network learning, as shown in Fig. 12.

4.1. Objective functions

The parameters of shafting alignment quality index include bearing support reaction force, shaft curve deflection and inter-shaft deflection. The parameters of shafting whirling vibration performance index include displacement amplitude, natural frequency and response amplitude of bearing force.

In order to simplify the calculation, the deflection of the shaft curve and the deflection angle between the segments are taken as computational constraints in the process of calculation. In order to avoid resonance and make the first-order natural frequency of the whirling vibration larger than the first critical speed of the shafting, the support reaction force at three bearings $F_i$ ($i=$1, 2, 3), the amplitude of three bearing monitoring points $L_i$ ($i=$1, 2, 3), the amplitude of the stern shaft monitoring point, the minimum safety factor of the shaft segment and first-order natural frequency $w_1$ in the multidisciplinary calculation model of the shafting are finally selected as the optimized objective variables, as shown in Table 2.

In order to improve the sensitivity of optimization parameters, the objective parameters are dimensionless, and the optimization rate index $\gamma$ is defined as follows:

where, $c_0$ is the initial value of the optimized parameters; $c_1$ is the value after the parameter optimization. the optimization rate is used to represent the optimization degree of the parameters in the table, so the dimensionless objective function can be further constructed:

where, ${\gamma }_j$ is the optimization variable numbered $j$ and ${\alpha }_j$is the optimization weight numbered $j$. According to the fuzzy comprehensive evaluation method based on Quality Function Deployment (QFD) [27-28], the weight matrix is given as follows:

Referring to the construction rules of fitness function in genetic algorithm, the fitness function is further obtained:

4.2. Optimization algorithm

According to the multi-disciplinary cross over model of shafting, fireworks algorithm was adopted as the optimization algorithm to accelerate the convergence of optimization and ensure the wholeness of optimization results.

The optimization mechanism of fireworks algorithm is similar to that of genetic algorithm, but the mechanism of controlling how much the sparks generated by fireworks within the explosion radius by explosion operator. The explosion radius makes the population with higher fitness produce more offspring in a smaller numerical range, and the population with lower fitness produce more dispersed and fewer offspring in a larger numerical range, which can better solve the optimal value search ability of the seven-dimensional design variables in this paper. The multi-disciplinary optimization process of Shafting based on fireworks algorithm is shown in Fig. 13.

Fig. 13Fireworks algorithm optimization process

Firstly, based on the principle of uniform randomness, the initial fireworks population was generated in the feasible region of each dimension design variable. After the preliminary calculation of the multi-disciplinary model of shafting, the fitness values of individuals in the initial population are calculated. The initial values of design variables and the feasible region are shown in Table 3.

The coding of each individual firework is set as follows in Fig. 14.

In Fig. 14, $X_i$ represents the vertical position coding section of bearing $i$, $Z_i$ represents the fulcrum position of bearing $i$, and $d_i$ represents the inner diameter of shaft $i$. For vertical position $X$, axial fulcrum position $Z$ and inner diameter $d$, the accuracy of 0.01mm is selected based on the actual technology and values. According to the coding rules, the coding section $X$ of each vertical position is 5 bits, the coding section $Z$ of axial position is 7 bits, and the coding section $d$ of inner radius is 4 bits.

Fig. 14The encoding of designs variable

Then, according to the explosion operator, the explosion radius $r_k$ and the number of sparks $N_k$ generated by the explosion of the $k$th firework are calculated. The sparks generated realize the migration and variation according to the migration rule $∆x$ and the mutation operator $∆E$. The formula of the explosion radius is expressed as follows:

where, ${\varphi }_{min}$ and ${\varphi }_{max}$ are the minimum and maximum fitness values in the current fireworks population, $\varepsilon$ is the minimum constant to avoid the denominator $r_k$ and $N_k$ being zero, $\varphi \left(k\right)$ is the fitness value of the $k$th fireworks, $K$ is the total number of fireworks, and the constants $R$ and $m$ are the magnitude of explosion radius and the maximum spark value, respectively.

Migration rules and mutation operators can be represented as follows:

where, $g\left(0,r_k\right)$ represents the random value of $r$ in [0, $r_k$], $\sigma$ represents the function value of Gaussian distribution, both the mean and variance are 1; $t$ represents the dimension serial number of an individual.

Finally, the migration and variation of cross-border design variables are mapped to the feasible region by mapping rules. The individuals with the maximum fitness value are retained according to the selection strategy, and the random individuals are retained by referring to spark distance, and finally, the next generation of fireworks is screened. Mapping rules and selection strategies based on spark distance can be expressed as follows:

In the above formula, $x_{max}^t$ and $x_{min}^t$ represent the upper and lower bounds of the individual on the dimension $t$ respectively, % represents modular operators, $d\left(x_k,x_j\right)$ represents the Euclidean distance between individuals, $D\left(x_k\right)$ represents the Euclidean distance between individuals, $P\left(x_k\right)$ represents the probability that each of the remaining individuals being kept randomly outside the best fitness individual, and $N$ represents the spark sets generated by explosion and mutation.

4.3. Optimization design results

Under the framework of multi-disciplinary optimization of shafting, the boundary conditions of optimization calculation should satisfy the following conditions:

1) For alignment calculation, each bearing load shall be within the upper limit [179.21 152.83 283.43] kN and the lower limit [11.43 7.85 6.28] kN;

2) The specific pressure of stern bearing shall not exceed 0.3MPa, and that of thrust bearing shall not exceed 0.25 MPa;

3) The safety factor of each section of the shafting should always be greater than 2.5;

4) The strength of each shaft segment shall meet the requirements, and the allowable stress ${\sigma }_{\mathrm{m}\mathrm{a}\mathrm{x}}\le \left[\sigma \right]$;

5) After optimization, the alignment technical parameters (e.g. deflection of shafting and rotation angle) and vibration technical parameters (e.g. anisotropic amplitude and critical speed), should meet the requirements in the specification [29-31].

The initial total number of fireworks is set as 2000 and the maximum spark is set as 1000 in the optimization algorithm. Through the joint simulation of MATLAB and finite element calculation, the internal diameter optimization value is 2.9 mm after iterative operation. The vertical displacement of rear stern bearing, front stern bearing and thrust bearing is [–8.89, 6.07, 0.72] mm, and the corresponding axial displacement is [–83.28, –44.6, –92, 55] respectively. Compared with the results of straight line alignment, the optimized results are shown in Fig. 14 and Table 4.

Fig. 15Whirling vibration frequency response curves of the bearing before and after MDO

a)$Y$-direction vibration of rear stern bearing before and after MDO

b)$Z$-direction vibration of rear stern bearing before and after MDO

c)$Y$-direction vibration of front stern bearing before and after MDO

d)$Y$-direction vibration of thrust bearing before and after MDO

From Fig. 15, it can be seen that the load of the rear stern bearing in the shafting has been improved after optimization with MDO method, compared with the initial position and size of the shafting. Within 0-50 Hz, the maximum vibration of the rear stern bearing and the front stern bearing has been improved. The maximum amplitude at the stern bearing and the front stern bearing decreases significantly, and the maximum amplitude of whirling vibration at the stern shaft decreased from 7.58 mm to 5.69 mm; the first 6 natural frequencies of the shafting slightly increase, the whirling vibration characteristics of the shafting are improved significantly.

Fig. 16Calibration results of shafting alignment after MDO

a) Deflection

b) Equivalent stress

By comparing Fig. 16 and Fig. 9, it can be seen that the maximum deflection of shafting increases from 0.4794 mm to 0.5759 mm, the maximum equivalent stress at the dangerous section decreases from 15.93 MPa to 15.06 MPa. It can be considered that the flexural characteristic of shafting is reduced, but the safety is increased. From Fig. 16 and Table 4, it can be seen that the load changes on each bearing are –7.72 kN, 22.16 kN, 9.01 kN, and the variations of load difference are 29.88 kN, 9.41 kN respectively. Although the loads of thrust bearing and front stern bearing are increased, the coincidence of the rear stern bearing is reduced, and the loading-distribution of the shafting is more uniform. It can be consider that the alignment characteristics of the shafting are improved.