## 2 Numerical Models  
### Governing Equation  
The water film pressure of a WLRSB considering surface roughness effects and 3D viscosity distribution can be determined by the 3D generalized average Reynolds equation, which can be expressed as:  
\[\frac{\partial}{R_{B}^{2}\partial\theta}\left(\phi_{\theta}G\frac{\partial p_{ h}}{\partial\theta}\right)+\frac{\partial}{\partial z}\left(\phi_{z}G_{R} \frac{\partial p_{h}}{\partial z}\right)=\omega\left(\phi_{c}\frac{\partial}{ \partial\theta}\left(h-\frac{I_{R}}{J_{R}}\right)+\frac{\sigma}{2}\frac{ \partial\phi_{s}}{\partial\theta}\right) \tag{1}\]  
where \(\phi_{\theta}\) and \(\phi_{z}\)[4] are the pressure flow factor in \(\theta\) and \(z\) directions, respectively; \(\phi_{c}\)[32] and \(\phi_{s}\)[4] are the contact factor and shear factor, respectively; \(p_{h}\) is the water film pressure, \(\rho_{W}\) is the density of water, and \(\eta\) is the viscosity of water. Furthermore, \(\omega\) is the speed of the journal, \(h\) is the average film thickness, and \(\sigma\) is the composite standard deviation of the rubber surface. In Equation (1), the parameters \(I_{R}\), \(J_{R}\), and \(G_{R}\) can be determined by:  
\[G_{R}=\int_{0}^{h}\frac{y}{\eta}\left(y-\frac{I_{R}}{J_{R}}\right)\mathrm{d}y,I_{R}=\int_{0}^{h}\frac{y}{\eta}\mathrm{d}y,J_{R}=\int_{0}^{h}\frac{1}{\eta} \mathrm{d}y \tag{2}\]  
where \(y\) indicates the radial coordinate along the cross-film direction. The profile of the WLRSB is shown in Figure 2. As shown, the lubrication clearance can be evaluated as the sum of the average geometric gap considering the deflection of the propeller shaft, \(h_{0}\), the groove depth, \(h_{g}\), and the elastic/thermal deformation of the shaft and the bush, \(\delta_{J}\) and \(\delta_{R}\). Consequently, the water film thickness, \(h(\theta,z,\Delta T)\), is given by the following:  
\[h(\theta,z,\Delta T)=h_{0}(\theta,z)+h_{g}(\theta,z)+\delta_{J}(\theta,\Delta T,p)+\delta_{R}(\theta,z,\Delta T,p) \tag{3}\]  
where \(z\) represents the axial direction shown in Figure 2, \(\Delta T\) is the temperature rise at the coordinate of \((\theta,z)\), and \(h_{0}\) can be calculated as:  
\[h_{0}=C(1+\varepsilon_{i}\cos(\theta-\varphi)) \tag{4}\]  
where \(\varepsilon_{i}\) is the eccentricity ratio of the \(i\)th cross-section of the bearings, and \(C\) is the radial clearance of the bearing. Figure 3 illustrates a schematic model of the deflected propeller shaft. In the force analysis of the propeller shaft, a simplification has been made so that the propeller shaft in the bearing hole is treated as a cantilever beam, which can be supported by two aspects. The first one is that the length-to-diameter ratio of the WLRSB is greatly larger than that of the general journal bearing, the second one is that the external load, i.e., the self-weight of the propeller, locates at the end face of the water-lubricated stern bearing. Additionally, a similar simplification was made in the study performed by [PERSON] et al. [8]. As shown, the bearing is discretized using \(n\) girds in the axial direction, the shaft is simplified as a cantilever beam, and the end face of the bearing that is close to the propeller is assumed to be fixed. The force acting on the \(i\)th section of the shaft is defined as \(F_{i}\), which is equal to the vector sum of \(F_{iz}\) in the \(z\) direction and \(F_{ix}\) in the \(x\) direction. It is noteworthy that both \(F_{iz}\) and \(F_{ix}\) comprise the hydrodynamic force component and the contact force component under FSH coupling conditions. The distances between the central point of \(i\)th section, \(O_{i}\), and point \(O^{\prime}\)\({}_{s}\) (which corresponds to the projection of the central point of the fixed journal face, \(O_{s}\), on the \(i\)th section A-A) are defined as \(d_{ix}\) in the \(x\) direction and \(d_{iy}\) in the \(y\) direction, as shown in Figure 3. According to the superposition principle, \(d_{ix}\) and \(d_{iy}\) can be calculated using the following formula:  
\[\left\{