On the Added Modal Coefficients of a Rotating Submerged Cylinder Induced by a Whirling Motion—Part 2:Numerical Investigation

Part 2 of this work presents a numerical methodology, validated using the experimental results presented in Part 1, to calculate the added modal coefficients of a submerged cylinder in water both when it oscillates and when it rotates with a whirling motion. The numerical methodology is based on computational fluid dynamic simulations that obtain the added modal forces on the cylinder when it is forced to vibrate with mode shapes calculated using acoustic-structural modal analysis. Then, these forces are processed with a curve-fitting algorithm to extract all the coefficients. Most numerical coefficients presented a close agreement with the corresponding experimental ones, although the added modal damping was overestimated. In general, the added modal mass was found to be independent of both the rotating speed and the whirling frequency except for low whirling frequencies when it increased. The added modal damping was found to depend on both parameters, and the rest of the coefficients were independent of the whirling frequency and only depended on the rotating speed. As a conclusion, this numerical approach has permitted the study of particular conditions that could not be experimentally tested and thus broadened the knowledge of the behavior of the added modal coefficients of rotating submerged cylinders.


Introduction
Most numerical studies devoted to calculating the added coefficients of non-rotating structures have been performed by means of coupled acoustic-structural simulations. This type of calculation consists of a two-way finite element simulation in which the solid mechanics are coupled with an acoustic domain representing the fluid. The acoustic elements model the fluid as irrotational and neglect the mean flow velocity. Consequently, the structural rotation cannot be considered in the analysis. Among all the added coefficients, these studies are limited to investigating the influence of the added mass on the modal response of non-rotating structures [1][2][3][4][5][6][7]. Escaler et al. [1] studied, both experimentally and numerically through coupled acoustic-structural simulations, the effects of water loading on the axisymmetric modes of vibration of a circular plate. A close agreement between experimental and numerical natural frequencies and mode shapes was found, with a frequency reduction ratio of around 64% due to the added mass effect. Although identical mode shapes were previously assumed for structures in air and fully submerged in water, measurable differences were experimentally and numerically observed in the radii of nodal circles between corresponding dry and wet modes. Similarly, Bossio et al. [2] also investigated the effects of the water surrounding a disk on its modal response by means of coupled acoustic-structural simulations. In addition to investigating the effects of the added mass, they studied how the acoustic natural frequencies of the fluid cavity alter the natural frequencies of the disk. The natural frequencies of the disk were found to decrease by up to 25% due to the proximity of an acoustic natural frequency. Moreover, the acoustic natural frequencies were observed to affect the disk natural frequencies only when both the acoustic and disk mode shapes presented the same number of nodal diameters. Liang et al. [3] investigated the added mass effects on the modal response of a Francis turbine runner using coupled acoustic-structural simulations. Compared with the corresponding experimental results, the numerical results showed strong agreement with a maximum deviation of around 3.5%. The frequency reduction ratio varied in the range from 10 to 39% depending on the mode shape. Huang et al. [4] also used acoustic-structural simulations to study how large-scale forms of attached cavitation, appearing on Francis runner blades, modified the added mass effects. It was concluded that the added mass effects induced on the entire runner were altered by the presence of the attached cavitation which was responsible for a decrease in the added mass effect and for an increase in the natural frequencies of the first modes. Moreover, any increase in cavitation cavity size was found to result in a decrement in the added mass effects and in an increment in the blade amplitude deformation below the cavity.
In order to calculate the added stiffness and damping of non-rotating structures, computational fluid dynamic (CFD) simulations which take into consideration the mean flow velocity are required. Liaghat et al. [8] performed two-way fluid-structure interaction (FSI) simulations to calculate the added damping using the exponential decay rate of the structure response to an impulse. Close agreements were found between experimental and numerical results with an average overestimation of the damping of around 12%. Monette et al. [9] implemented a theoretical model in an in-house finite element (FE) code based on kinetic energy transfer between the fluid flow and the moving structure to calculate the added damping. A strong correlation was found between experimental and numerical results. Both researchers relied on a two-way coupling of the fluid flow and the structural dynamics, but this implementation is complex and requires a high computational effort. To overcome this problem, Gauthier et al. [10] proposed a method to calculate the added mass, damping, and stiffness based on a prescribed motion of the structural boundary in Unsteady Reynolds-Averaged Navier-Stokes (URANS) simulations. This method does not require two-way simulations; nonetheless, the natural frequency of the coupled system must be known beforehand. The natural frequency and added mass were calculated using an acoustic-structural modal analysis prior to performing: (i) a URANS simulation to calculate the added damping, and (ii) a Reynolds-Averaged Navier-Stokes (RANS) simulation to calculate the added stiffness. Roig et al. [11] used a similar numerical methodology to investigate the influence of the wake cavitation on the dynamic response of hydraulic profiles under lock-in conditions. The numerical results correlated well with the corresponding experimental results, which also validated the methodology.
A non-rotating structure has standing mode shapes whilst a rotating structure has whirling mode shapes. Consequently, the added coefficients of the submerged components of rotating machines with a whirling motion, such as runners, should be considered in their dynamic analysis. The few researchers who have approached this topic have only studied the influence of the runner's added mass on the dynamic response of rotors without considering other potential added coefficients. Moreover, the added mass exerted on the runner has been calculated by means of acoustic-structural simulations and therefore without considering the influence of the rotation or whirling motion. Gustavsson et al. [12] calculated the added mass and moments of inertia of turbine runners using acousticstructural simulations. Their values were then introduced in a rotor dynamic simulation. The added mass of the runner was found to be influenced by its geometry, blade angle, and the shape of the draft tube walls. It was observed that both the polar moment of inertia and the radial mass increased by about 300% and 60%, respectively. Similarly, Keto-Tokio et al. [13] also calculated the added mass of a turbine runner through acoustic-structural simulations and included it in a rotor dynamic analysis. The lowest natural frequencies of the turbine shaft train were found to decrease by 25%.
To the author's knowledge, only numerical methodologies to calculate the added coefficients of structures, such as hydraulic machines, without rotation have been presented. Consequently, it is necessary to develop and validate numerical methodologies to calculate and investigate the added coefficients of submerged bodies induced by the rotation and whirling motion.
This article presents a numerical methodology, validated using the experimental results of Part 1 [14], to calculate the added modal coefficients of a submerged cylinder both when it only oscillates and when it rotates and presents a whirling motion. The test rig and the experimental methodology used to measure and calculate the added modal coefficients are presented in Part 1 [14]. The numerical methodology breaks down the problem leading to a deeper understanding of the physics involved. Furthermore, the presented numerical results provide an insight into how and why added modal coefficients are influenced by the rotating speed and the whirling frequency.

Numerical Methodology
The following presents a numerical methodology to calculate: (i) the added modal mass, M f , and the added modal damping, C f , of a standing mode shape when the cylinder does not rotate, and (ii) the added modal mass, M f , the added modal damping, C f , the added modal stiffness, K f , and the two new added modal coefficients proposed in Part 1 [14], A f and B f , of a whirling mode shape when the cylinder rotates. One advantage of the numerical model compared to the experiments carried out and presented in Part 1 [14] is that the mode shapes can be forced to oscillate at any selected frequency, without being restricted to the cylinder natural frequencies.
The first step of this method is to map the normalized structural mode shape of interest on the fluid-structure interface of a URANS simulation. Subsequently, the mapped mode shape is forced to only oscillate, under non-rotating conditions, or to oscillate and rotate, under rotating conditions. Finally, the added modal coefficients can be obtained using the calculated added modal forces exerted by the fluid on the cylinder.
Assuming that the mean flow velocity does not alter the cylinder mode shapes, they can be calculated through an acoustic-structural modal analysis. Figure 1 shows the first two mode shapes of the cylinder in water before unity normalization, Φ f 1 and Φ f 2 , which have the same natural frequency and are identical but perpendicular to each other.  To the author's knowledge, only numerical methodologies to calculate the ad  efficients of structures, such as hydraulic machines, without rotation have been pr  Consequently, it is necessary to develop and validate numerical methodologies t  late and investigate the added coefficients of submerged bodies induced by the  and whirling motion. This article presents a numerical methodology, validated using the experim sults of Part 1 [14], to calculate the added modal coefficients of a submerged cylin when it only oscillates and when it rotates and presents a whirling motion. The and the experimental methodology used to measure and calculate the added mod ficients are presented in Part 1 [14]. The numerical methodology breaks down the leading to a deeper understanding of the physics involved. Furthermore, the pr numerical results provide an insight into how and why added modal coefficient fluenced by the rotating speed and the whirling frequency.

Numerical Methodology
The following presents a numerical methodology to calculate: (i) the added mass, , and the added modal damping, , of a standing mode shape when the does not rotate, and (ii) the added modal mass, , the added modal damping added modal stiffness, , and the two new added modal coefficients proposed [14], and , of a whirling mode shape when the cylinder rotates. One advantag numerical model compared to the experiments carried out and presented in Part that the mode shapes can be forced to oscillate at any selected frequency, witho restricted to the cylinder natural frequencies.
The first step of this method is to map the normalized structural mode shap terest on the fluid-structure interface of a URANS simulation. Subsequently, the mode shape is forced to only oscillate, under non-rotating conditions, or to oscil rotate, under rotating conditions. Finally, the added modal coefficients can be o using the calculated added modal forces exerted by the fluid on the cylinder.
Assuming that the mean flow velocity does not alter the cylinder mode shap can be calculated through an acoustic-structural modal analysis. Figure 1 shows two mode shapes of the cylinder in water before unity normalization, and have the same natural frequency and are identical but perpendicular to each othe Under non-rotating conditions, was normalized, imported, and mapped fluid-structure interface of a URANS simulation. The mapped normalized mod Under non-rotating conditions, Φ f 1 was normalized, imported, and mapped on the fluid-structure interface of a URANS simulation. The mapped normalized mode shape was forced to oscillate at a selected frequency, f , and with a selected modal amplitude, h, which scaled the normalized mode shape. The added modal force, F f q 1 q 1 , induced as a reaction to the modal acceleration ( .. q 1 ), velocity ( . q 1 ) and displacement (q 1 ) of the cylinder: .. .
in the direction associated with Φ f 1 was then calculated. Once the solution was stable, eight periods of F f q 1 q 1 were processed by a curve-fitting algorithm which estimated F f q 1 q 1 through a harmonic signal: by a least square method (see Figure 2). As shown in Figure 2, it must be noted that under non-rotating conditions the fluid only exerts a force in the direction associated with Φ f 1 , whereas the force in the direction associated with Φ f 2 , F f q 2 q 1 , is zero. was forced to oscillate at a selected frequency, , and with a selected modal amplitude, ℎ, which scaled the normalized mode shape. The added modal force, , induced as a reaction to the modal acceleration ( ), velocity ( ) and displacement ( ) of the cylinder: in the direction associated with was then calculated. Once the solution was stable, eight periods of were processed by a curve-fitting algorithm which estimated through a harmonic signal: by a least square method (see Figure 2). As shown in Figure 2, it must be noted that under non-rotating conditions the fluid only exerts a force in the direction associated with , whereas the force in the direction associated with , , is zero. Assuming that under non-rotating conditions a static displacement, , of the submerged cylinder does not induce per se added forces, can be modeled as a function of and as well as their corresponding added modal coefficients, as shown in Equation (5). and were then calculated from a direct comparison of and , which were computed using the numerical model and the curve-fitting algorithm, with Equation (5) (see Equations (6) and (7)).  Assuming that under non-rotating conditions a static displacement, q 1 , of the submerged cylinder does not induce per se added forces, F f q 1 q 1 can be modeled as a function of .. q 1 and . q 1 as well as their corresponding added modal coefficients, as shown in Equation (5). M f and C f were then calculated from a direct comparison of α and β, which were computed using the numerical model and the curve-fitting algorithm, with Equation (5) (see Equations (6) and (7)).
Under rotating conditions, the whirling mode shapes can be understood as the combination of Φ f 1 and Φ f 2 , oscillating with the same h and at the same f but at 90 • out of phase relative to each other. Thus, having ..  q 2 , and q 2 associated with Φ f 2 must be defined as: .. .
In this case, there are two forces: (i) F f q 1 in the direction associated with Φ f 1 induced by ..
. q 2 , q 1 , and q 2 as well as their corresponding added modal coefficients, and (ii) F f q 2 in the direction associated with Φ f 2 induced by ..
. q 1 , q 2 , and q 1 as well as their corresponding added modal coefficients (see Equations (11) and (12)): Nonetheless, the whirling motion was simplified by only considering the oscillation of Φ f 1 and the rotation rather than considering the coupling of both Φ f 1 and Φ f 2 , which would reproduce the full whirling motion. The use of only one mode shape simplified the numerical setup, reduced the number of simulations to calculate all the added modal coefficients, and minimized the required computation resources, whereas the final results appeared to be accurate. Considering this simplification, F f q 1 q 1 in the direction associated with Φ f 1 is only induced by .. q 1 , . q 1 , and q 1 as well as their corresponding added modal coefficients, and F f q 2 q 1 in the direction associated with Φ f 2 is induced by . q 1 and q 1 as well as their corresponding added modal coefficients (see Equations (13) and (14)): As shown qualitatively by the black and red arrows in Figure 3, the sine amplitude of F f q 1 is approximately equal to the sine amplitude of F f q 1 q 1 plus the absolute value of the cosine amplitude of F f q 2 q 1 . Using the curve-fitting algorithm, the sine coefficient of F f q 1 can be calculated and is 1.23 N. Similarly, the sine coefficient of F f q 1 q 1 and the absolute value of the cosine coefficient of F f q 2 q 1 can also be calculated and are 0.7563 and |−0.4801| N, respectively, which when added give a value of about 1.2364 N showing a deviation of 0.5% relative to the actual amplitude, calculated using the full whirling motion. Thus, it is verified that the simplification is approximately valid. Both simulated and were approximated by the curve-fitting algorithm and their coefficients compared with Equations (13) and (14), as shown in Equations (15) and (16) for the case of : and in Equations (17) and (18) for the case of : and were calculated directly from Equations (16) and (18), respectively. In order to separate the contributions of and in Equation (15), a new simulation was performed which consisted of forcing the fluid-structure interface to rotate and deform statically according to scaled by ℎ. Under these circumstances, in the direction associated with is only induced by as well as its corresponding added modal coefficient, and in the direction associated with is only induced by as well as its corresponding added modal coefficient, as shown in Equations (19) and (20) and Figure  4: ℎ, Both simulated F f q 1 q 1 and F f q 2 q 1 were approximated by the curve-fitting algorithm and their coefficients compared with Equations (13) and (14), as shown in Equations (15) and (16) for the case of F f q 1 q 1 : (16) and in Equations (17) and (18) for the case of F f q 2 q 1 : C f and A f were calculated directly from Equations (16) and (18), respectively. In order to separate the contributions of M f and K f in Equation (15), a new simulation was performed which consisted of forcing the fluid-structure interface to rotate and deform statically according to Φ f 1 scaled by h. Under these circumstances, F f q 1 q 1 in the direction associated with Φ f 1 is only induced by q 1 as well as its corresponding added modal coefficient, and F f q 2 q 1 in the direction associated with Φ f 2 is only induced by q 1 as well as its corresponding added modal coefficient, as shown in Equations (19) and (20) and Figure 4: K f and B f were subsequently computed substituting the numerically calculated values of F f q 1 q 1 and F f q 2 q 1 in Equations (19) and (20). Finally, using K f and Equation (15), M f was computed.
It must be noted that under rotating conditions the numerical forces present a ripple (see Figure 3) which is induced by the mesh when it rotates and deforms, but is not a physical phenomenon. and were subsequently computed substituting the numerically calculated values of and in Equations (19) and (20). Finally, using and Equation (15), was computed.
It must be noted that under rotating conditions the numerical forces present a ripple (see Figure 3) which is induced by the mesh when it rotates and deforms, but is not a physical phenomenon.

Acoustic-Structural Modal Analysis
Coupled acoustic-structural simulations have been widely used to calculate natural frequencies and mode shapes of submerged bodies showing close agreement with corresponding experimental results [1][2][3][4][5][6][7]. This method usually consists of coupling structural elastic elements with potential flow elements based on the Laplace equation or with acoustic elements based on the Helmholtz equation. In the present work, the coupled acousticstructural simulations were conducted with the commercial solver code Ansys which models the fluid by acoustic elements. Based on the theory presented by [15], the continuity equations and the momentum equations can be modified assuming that: (i) the fluid is compressible, (ii) the fluid is irrotational, (iii) there is no body force, (iv) the pressure disturbance of the fluid is small, (v) there is no mean flow of the fluid, and (vi) the gas is ideal, adiabatic and reversible, resulting in the following linearized continuity and momentum equations: where and are the acoustic velocity and pressure, respectively, is a mass source, is the speed of sound in the fluid, is the mean fluid density, is the dynamic viscosity of the fluid, and is time. Finally, the acoustic wave equation can be given by:

Acoustic-Structural Modal Analysis
Coupled acoustic-structural simulations have been widely used to calculate natural frequencies and mode shapes of submerged bodies showing close agreement with corresponding experimental results [1][2][3][4][5][6][7]. This method usually consists of coupling structural elastic elements with potential flow elements based on the Laplace equation or with acoustic elements based on the Helmholtz equation. In the present work, the coupled acoustic-structural simulations were conducted with the commercial solver code Ansys which models the fluid by acoustic elements. Based on the theory presented by [15], the continuity equations and the momentum equations can be modified assuming that: (i) the fluid is compressible, (ii) the fluid is irrotational, (iii) there is no body force, (iv) the pressure disturbance of the fluid is small, (v) there is no mean flow of the fluid, and (vi) the gas is ideal, adiabatic and reversible, resulting in the following linearized continuity and momentum equations: where v a and p a are the acoustic velocity and pressure, respectively, Q is a mass source, c is the speed of sound in the fluid, ρ f is the mean fluid density, µ is the dynamic viscosity of the fluid, and t is time. Finally, the acoustic wave equation can be given by: The finite element formulation is obtained by testing Equation (23) using the Galerkin procedure [16] and combining it with the equation of the normal velocity, v n,F , on the boundary of the acoustic domain: where u F is the displacement of fluid particles andn is the normal vector. Describing the pressure and displacement components using shape functions, and assuming that there are no fluid forces, the discretized acoustic wave presented in Equation (23) can be expressed in matrix notation as: which can be coupled with the structural dynamic equation of a body submerged in an acoustic domain, which takes the following form when assuming that there are no structural forces:  (27), the natural frequencies and mode shapes of bodies submerged in a fluid taking into consideration the added mass effects can be computed. Among all the parts of the test rig, only the slim shaft and the cylinder were included in the acoustic-structural model since the mode shapes of interest only affected this region, as shown in Figure 1 where it is observed that the top of the slim shaft presents zero displacement. The material of the shaft and the cylinder was considered a standard stainless steel. The acoustic domain was modeled as water with a density and speed of sound of 1000 kg/m 3 and 1430 m/s, respectively. Showing the boundary conditions in Figure 5, it can be seen that: (i) the top of the shaft was considered as clamped, (ii) the top and bottom covers as well as the lateral walls of the cylindrical water tank were considered rigid walls, and (iii) the nodes of the solid elements in contact with the fluid were considered as part of a fluid-solid interface.
The finite element formulation is obtained by testing Equation (23) using the Galerkin procedure [16] and combining it with the equation of the normal velocity, , , on the boundary of the acoustic domain: where is the displacement of fluid particles and is the normal vector. Describing the pressure and displacement components using shape functions, and assuming that there are no fluid forces, the discretized acoustic wave presented in Equation (23) can be expressed in matrix notation as: which can be coupled with the structural dynamic equation of a body submerged in an acoustic domain, which takes the following form when assuming that there are no structural forces: , (27) where and are the acoustic and structural mass matrices, and are the acoustic and structural damping matrices, and are the acoustic and structural stiffness matrices, is a coupling matrix, and is the structural displacement which is equal to at the fluid-structure interface. Solving Equation (26) coupled with Equation (27), the natural frequencies and mode shapes of bodies submerged in a fluid taking into consideration the added mass effects can be computed.
Among all the parts of the test rig, only the slim shaft and the cylinder were included in the acoustic-structural model since the mode shapes of interest only affected this region, as shown in Figure 1 where it is observed that the top of the slim shaft presents zero displacement. The material of the shaft and the cylinder was considered a standard stainless steel. The acoustic domain was modeled as water with a density and speed of sound of 1000 kg/m 3 and 1430 m/s, respectively. Showing the boundary conditions in Figure 5, it can be seen that: (i) the top of the shaft was considered as clamped, (ii) the top and bottom covers as well as the lateral walls of the cylindrical water tank were considered rigid walls, and (iii) the nodes of the solid elements in contact with the fluid were considered as part of a fluid-solid interface. The discretization of the fluid and structural domains was performed using tetrahedral elements. A mesh refinement process was carried out to determine the optimal element size that gives an accurate solution with the lowest computational cost. The value of the natural frequency was considered as the verification variable and it was plotted as a function of the number of mesh nodes (see Figure 6). The final optimal mesh had 962,609 nodes. The discretization of the fluid and structural domains was performed using tetrahedral elements. A mesh refinement process was carried out to determine the optimal element size that gives an accurate solution with the lowest computational cost. The value of the natural frequency was considered as the verification variable and it was plotted as a function of the number of mesh nodes (see Figure 6). The final optimal mesh had 962,609 nodes. Figure 6. Mesh refinement study for the acoustic-structural modal analysis.

Computational Fluid Dynamic Analysis
The CFD simulations were conducted with the commercial solver Ansys CFX. The fluid was described using the continuity and momentum equations which were averaged when solving a turbulent flow. The turbulent flow variables are separated into average and time-varying components based on the theory presented by [17]. If the fluid is assumed as incompressible and the body forces are neglected, the following conservation equations are obtained: where and are the average and varying components of the velocity, respectively, is the average pressure, is the fluid kinematic viscosity, and represents the Reynolds stress tensor.
The unknown Reynolds stresses are obtained from the Boussinesq hypothesis: and using the shear stress transport (SST) model which combines the k-ω model near the walls and the k-ε model away from the walls. The SST and other similar turbulence models have already been used in the research of whirling flows with excellent results [18][19][20]. The turbulent viscosity, , can be calculated as: Figure 6. Mesh refinement study for the acoustic-structural modal analysis.

Computational Fluid Dynamic Analysis
The CFD simulations were conducted with the commercial solver Ansys CFX. The fluid was described using the continuity and momentum equations which were averaged when solving a turbulent flow. The turbulent flow variables are separated into average and time-varying components based on the theory presented by [17]. If the fluid is assumed as incompressible and the body forces are neglected, the following conservation equations are obtained: where V and v are the average and varying components of the velocity, respectively, P is the average pressure, ν is the fluid kinematic viscosity, and ρ f v i v j represents the Reynolds stress tensor. The unknown Reynolds stresses are obtained from the Boussinesq hypothesis: and using the shear stress transport (SST) model which combines the k-ω model near the walls and the k-ε model away from the walls. The SST and other similar turbulence models have already been used in the research of whirling flows with excellent results [18][19][20]. The turbulent viscosity, ν t , can be calculated as: where k and ω are obtained by solving the following transport equations: the blending function F 1 is expressed as: and CD kω as: where y is the distance to the nearest wall, S is the invariant strain rate, and F 2 is defined as: The turbulence energy is limited through: where P k is defined as: All constants are computed via: and have the following empirical values: (i) σ k1 = 0.85, (ii) σ k2 = 1.00, (iii) σ ω1 = 0.50, (iv) σ ω2 = 0.856, (v) β * = 0.09, (vi) β 1 = 0.075, (vii) β 1 = 0.0828, (viii) α 1 = 5/9, (ix) α 2 = 0.44. The flow was modeled as laminar for the non-rotating cases since the lowest and highest Reynolds numbers were 38.5 and 2705, respectively. The Reynolds number was defined based on the diameter of the cylinder, 0.164 m, and the average between the maximum velocities presented by the cylinder along its span. Nevertheless, when the cylinder rotates, the flow has been modeled as turbulent since the lowest and highest Reynolds numbers among the different cases studied were 3864 and 48,302, respectively. Under rotating conditions, the Reynolds number was defined based on the gap between the cylinder and the wall, 0.015 m, and the maximum flow velocity in the gap, induced only by the cylinder rotation.
In this sense, a high-resolution turbulence model was used to solve the flow when the cylinder rotates which resolves smaller turbulent scales than those solved by a first-order turbulence model because the former model uses the high-resolution advection scheme and the second-order backward Euler transient scheme.
The CFD domain was modeled using the water material defined in Section 2.1.1. Showing the boundary conditions in Figure 7, it can be seen that: (i) the top and bottom covers as well as the lateral walls of the cylindrical water tank were defined as stationary no-slip walls, (ii) the fluid-structure interface was also defined as a no-slip wall but it was forced to oscillate periodically according to Φ f1 , and (iii) under rotating conditions, the rotating subdomain was also forced to rotate whilst a general connection interface (GCI) was used to join both subdomains. forced to oscillate periodically according to , and (iii) under rotating conditions, the rotating subdomain was also forced to rotate whilst a general connection interface (GCI) was used to join both subdomains. The validation of the grid independence was carried out under non-rotating conditions considering as a reference the values of and . Figure 8 shows and as a function of the number of nodes when the mesh resolution was increased in the axial, circumferential, and radial directions.
shows approximately the same value with all the mesh sizes. Nonetheless, increases and decreases when the number of nodes in radial directions as well as in axial and circumferential directions is increased.  The validation of the grid independence was carried out under non-rotating conditions considering as a reference the values of M f and C f . Figure 8 shows M f and C f as a function of the number of nodes when the mesh resolution was increased in the axial, circumferential, and radial directions. M f shows approximately the same value with all the mesh sizes.
Nonetheless, C f increases and decreases when the number of nodes in radial directions as well as in axial and circumferential directions is increased. forced to oscillate periodically according to , and (iii) under rotating conditions, the rotating subdomain was also forced to rotate whilst a general connection interface (GCI) was used to join both subdomains. The validation of the grid independence was carried out under non-rotating conditions considering as a reference the values of and . Figure 8 shows and as a function of the number of nodes when the mesh resolution was increased in the axial, circumferential, and radial directions.
shows approximately the same value with all the mesh sizes. Nonetheless, increases and decreases when the number of nodes in radial directions as well as in axial and circumferential directions is increased.   Figure 9 presents the results of a sensitivity study of the time step, t s , to obtain an independent model. M f and C f are plotted as a function of the number of time steps per vibration period, T, through T/t s . It is observed that M f presents the same value for all t s ; however, C f decreases when t s is decreased.
J. Mar. Sci. Eng. 2023, 11, x FOR PEER REVIEW 12 of 23 Figure 9 presents the results of a sensitivity study of the time step, ts, to obtain an independent model. and are plotted as a function of the number of time steps per vibration period, T, through T/ts. It is observed that presents the same value for all ts; however, decreases when ts is decreased. It can be observed that for the smallest time step considered corresponding to /1000, a numerical result of 5.6 Ns/m was obtained which was quite different from the experimental value of 4.23 Ns/m. Due to the excessive computational time required to achieve a valid solution for such short ts, it was decided to select a mesh with 73,748 nodes for a /100 as it was sufficient to provide a mesh independent value of but it was assumed to work with numerical values of significantly higher than the experimental ones. Figure 10a shows the selected structured mesh built using hexahedral elements with a detail of the mesh around the cylinder, which gives approximately 16 mesh layers in the boundary layer and a Y+ lower than 1 as shown in Figure 10b.  It can be observed that for the smallest time step considered corresponding to t s = T/ 1000, a numerical result of C f = 5.6Ns/m was obtained which was quite different from the experimental value of 4.23 Ns/m. Due to the excessive computational time required to achieve a valid solution for such short t s , it was decided to select a mesh with 73,748 nodes for a t s = T/100 as it was sufficient to provide a mesh independent value of M f but it was assumed to work with numerical values of C f significantly higher than the experimental ones. Figure 10a shows the selected structured mesh built using hexahedral elements with a detail of the mesh around the cylinder, which gives approximately 16 mesh layers in the boundary layer and a Y+ lower than 1 as shown in Figure 10b.  Figure 9 presents the results of a sensitivity study of the time step, ts, to obtain an independent model. and are plotted as a function of the number of time steps per vibration period, T, through T/ts. It is observed that presents the same value for all ts; however, decreases when ts is decreased. It can be observed that for the smallest time step considered corresponding to /1000, a numerical result of 5.6 Ns/m was obtained which was quite different from the experimental value of 4.23 Ns/m. Due to the excessive computational time required to achieve a valid solution for such short ts, it was decided to select a mesh with 73,748 nodes for a /100 as it was sufficient to provide a mesh independent value of but it was assumed to work with numerical values of significantly higher than the experimental ones. Figure 10a shows the selected structured mesh built using hexahedral elements with a detail of the mesh around the cylinder, which gives approximately 16 mesh layers in the boundary layer and a Y+ lower than 1 as shown in Figure 10b.

Case Studies
Multiple simulations have been performed with always the same h and multiple combinations of rotating speeds, Ω, and f , which have been summarized in Table 1.

Numerical Methodology Validation
M f , C f , A f and K f were calculated numerically for different Ω and f to be compared with the experiments performed in Part 1 [14] (see Table 1). Figure 11 shows a comparison between the numerical (Num) and experimental (Exp) added modal coefficients. The numerical M f calculated when the cylinder does not rotate coincides with the experimental one, presenting a deviation of about 1.2%. The numerical C f always presents higher values than the experimental ones, but its trend is well captured at low Ω. The average deviation is 142.5% and may be partially induced by an additional damping added by the mesh when it deforms and by the use of a not sufficiently short t s . In previous studies, a numerical damping higher than the corresponding experimental one was also obtained [8,21]. The numerical and experimental A f and K f present a strong agreement with average deviations around 10% and 18%, respectively. In summary, it is proved that the numerical model provides good accuracy for most of the added modal coefficients with the exception of C f .

Case Studies
Multiple simulations have been performed with always the same ℎ and multiple combinations of rotating speeds, , and , which have been summarized in Table 1.

Numerical Methodology Validation
, , and were calculated numerically for different and to be compared with the experiments performed in Part 1 [14] (see Table 1). Figure 11 shows a comparison between the numerical (Num) and experimental (Exp) added modal coefficients. The numerical calculated when the cylinder does not rotate coincides with the experimental one, presenting a deviation of about 1.2%. The numerical always presents higher values than the experimental ones, but its trend is well captured at low . The average deviation is 142.5% and may be partially induced by an additional damping added by the mesh when it deforms and by the use of a not sufficiently short . In previous studies, a numerical damping higher than the corresponding experimental one was also obtained [8,21]. The numerical and experimental and present a strong agreement with average deviations around 10% and 18%, respectively. In summary, it is proved that the numerical model provides good accuracy for most of the added modal coefficients with the exception of .  Figure 12 shows the evolution of and as a function of . presents a constant value when is higher than 1 Hz, which agrees with previous studies that concluded that was independent of the flow conditions and [10,[22][23][24]. However, for lower than 1 Hz, the increases when is decreased, as was also observed in a previous investigation with Kaplan turbines [25].

Added Modal Coefficients under Non-Rotating Conditions
presents an approximately linear increase when is increased. Nonetheless, for lower than 4.5 Hz, deviates from the linear trend, presenting higher values than those expected. It can be seen that evolves similarly to a drag force, which increases quadratically with the flow velocity, when assuming that the flow velocity is equivalent to . Consequently, can be understood as a drag force type.  Table 2 shows the percentage of and induced by the pressure (PFs) and viscous (VFs) forces as a function of . In both cases, the highest contribution of the VFs to and occurs for the lowest , representing a contribution of 1.09% and 6.94%, respectively. At higher values of , its contribution decreases. For instance, at 7 Hz, VFs only contribute 0.14% to and 2.92% to . PFs are the main contributors to and due to the cylindrical geometry, especially its large cross-section area. Finally, it can be seen that the contribution of VFs to is always higher than to .  Figure 12 shows the evolution of M f and C f as a function of f . M f presents a constant value when f is higher than 1 Hz, which agrees with previous studies that concluded that M f was independent of the flow conditions and f [10,[22][23][24]. However, for f lower than 1 Hz, the M f increases when f is decreased, as was also observed in a previous investigation with Kaplan turbines [25].  Figure 12 shows the evolution of and as a function of . presents a constant value when is higher than 1 Hz, which agrees with previous studies that concluded that was independent of the flow conditions and [10,[22][23][24]. However, for lower than 1 Hz, the increases when is decreased, as was also observed in a previous investigation with Kaplan turbines [25].

Added Modal Coefficients under Non-Rotating Conditions
presents an approximately linear increase when is increased. Nonetheless, for lower than 4.5 Hz, deviates from the linear trend, presenting higher values than those expected. It can be seen that evolves similarly to a drag force, which increases quadratically with the flow velocity, when assuming that the flow velocity is equivalent to . Consequently, can be understood as a drag force type.  Table 2 shows the percentage of and induced by the pressure (PFs) and viscous (VFs) forces as a function of . In both cases, the highest contribution of the VFs to and occurs for the lowest , representing a contribution of 1.09% and 6.94%, respectively. At higher values of , its contribution decreases. For instance, at 7 Hz, VFs only contribute 0.14% to and 2.92% to . PFs are the main contributors to and due to the cylindrical geometry, especially its large cross-section area. Finally, it can be seen that the contribution of VFs to is always higher than to . C f presents an approximately linear increase when f is increased. Nonetheless, for f lower than 4.5 Hz, C f deviates from the linear trend, presenting higher values than those expected. It can be seen that C f . q 1 evolves similarly to a drag force, which increases quadratically with the flow velocity, when assuming that the flow velocity is equivalent to . q 1 . Consequently, C f . q 1 can be understood as a drag force type. Table 2 shows the percentage of M f and C f induced by the pressure (PFs) and viscous (VFs) forces as a function of f . In both cases, the highest contribution of the VFs to M f and C f occurs for the lowest f , representing a contribution of 1.09% and 6.94%, respectively. At higher values of f , its contribution decreases. For instance, at 7 Hz, VFs only contribute 0.14% to M f and 2.92% to C f . PFs are the main contributors to M f and C f due to the cylindrical geometry, especially its large cross-section area. Finally, it can be seen that the contribution of VFs to C f is always higher than to M f .
where A is the cylinder cross-section area. It can be seen that M f * and C f * tend to a constant value when f is higher than 1 and 4.5 Hz, respectively, and that they increase when f is decreased at lower values.
where is the cylinder cross-section area. It can be seen that * and * tend to a constant value when is higher than 1 and 4.5 Hz, respectively, and that they increase when is decreased at lower values.
(a) (b)  Figure 14a shows an isometric view of the cylinder with a green line over the perimeter located at the midplane. The red arrow indicates the direction and sense of the vibration at time instants , corresponding to the maximum amplitude of the oscillation, and 0 , corresponding to the zero oscillating amplitude. It must be noted that at , the cylinder presents the maximum and but a 0 m/s. On the other hand, at 0 , the cylinder presents the maximum but a 0 m/s 2 and a 0 m.
Consequently, at and 0 the added modal forces are induced exclusively by and , respectively. Figure 14b shows the locations on the green perimeter corresponding to the axis of Figures 15 and 16, with 0 rad corresponding to the front of the cylinder and rad to the back of the cylinder.  Figure 14a shows an isometric view of the cylinder with a green line over the perimeter located at the midplane. The red arrow indicates the direction and sense of the vibration at time instants t = 3 4 T, corresponding to the maximum amplitude of the oscillation, and t = 0T, corresponding to the zero oscillating amplitude. It must be noted that at t = 3 4 T, the cylinder presents the maximum ..    (a) (b) Figure 14. Isometric view of the cylinder with a green line over the perimeter located at the midplane (a) and reference locations on the green perimeter (b). Figure 15 shows the normalized pressure plotted along the green perimeter at induced by , * , and at 0 induced by , * , as: At 0 rad, both * and * present the highest value as shown in Figure 15a,b. * and * decrease and become negative when moving along the perimeter in the counter-clockwise direction. At rad , * and * show the lowest pressure and then increase progressively when moving towards to the front of the cylinder. It must be noted that lower results in higher * and * and that at higher values of , they tend to present the same amplitudes regardless of .
(a) (b) Figure 15. Normalized pressure plotted along the green perimeter at (a) and at 0 (b). Figure 16 shows the normalized velocity gradient in Y direction, * , plotted along the green perimeter at and consequently related exclusively with as: where is the velocity gradient in direction Y. As shown in Figure 16, at 0 rad, * 0 and then decreases and becomes negative moving over the perimeter in the counter-clockwise direction, reaching the lowest value at rad . Subsequently, * increases progressively up to 0 at rad , then continues increasing and presents the highest value at rad. Finally, * decreases moving towards to the front of the cylinder. It must be noted that lower values of result in lower | * | and that at higher values of , it tends to present the same amplitudes regardless of . Based on Figures 12b, 13b, and 15b, it can be concluded that loses its linearity at low values of , and thus when there is the highest influence of the VFs on (see Table  2). Consequently, the increase in at low values of may be partially induced by this  Figure 15 shows the normalized pressure plotted along the green perimeter at t = 3 4 T induced by M f , P M f * , and at t = 0T induced by C f , P C f * , as: At θ = 0 rad, both P M f * and P C f * present the highest value as shown in Figure 15a,b. P M f * and P C f * decrease and become negative when moving along the perimeter in the counter-clockwise direction. At θ = π rad, P M f * and P C f * show the lowest pressure and then increase progressively when moving towards to the front of the cylinder. It must be noted that lower f results in higher P M f * and P C f * and that at higher values of f , they tend to present the same amplitudes regardless of f . Figure 16 shows the normalized velocity gradient in Y direction, V Y * , plotted along the green perimeter at t = 3 4 T and consequently related exclusively with M f as: where V Y is the velocity gradient in direction Y. As shown in Figure 16, at θ = 0 rad, V Y * = 0 and then decreases and becomes negative moving over the perimeter in the counter-clockwise direction, reaching the lowest value at θ = π 2 rad. Subsequently, V Y * increases progressively up to 0 at θ = π rad, then continues increasing and presents the highest value at θ = 3π 2 rad. Finally, V Y * decreases moving towards to the front of the cylinder. It must be noted that lower values of f result in lower |V Y * | and that at higher values of f , it tends to present the same amplitudes regardless of f . Based on Figure 12b, Figure 13b, and Figure 15b, it can be concluded that C f loses its linearity at low values of f , and thus when there is the highest influence of the VFs on C f (see Table 2). Consequently, the increase in C f at low values of f may be partially induced by this increase in the influence of VFs.
Based on Figure 12a, Figure 13a, and Figure 15a, it can be concluded that M f increases at low f when the V Y * presents the lowest values (see Figure 16). Low values of V Y * suggest that the velocity over all the length of the gap between the cylinder and the tank wall is approximately equal, indicating that the acceleration of the cylinder accelerates all the fluid of the gap rather than only the fluid close to the cylinder as at high values of V Y * . It must result in a higher amount of water mass accelerated and thus in higher M f . This may be partially induced by the increase in the influence of VFs at low f (see Table 2). Figures 17 and 18 show the distributions of the pressure and velocity, respectively, on the midplane of the cylindrical water tank when Ω = 1.25 and 6.25 Hz and f = 0 Hz. Under these circumstances, the cylinder is displaced, according to the Φ f1 scaled by h, and it does not oscillate but only rotates. The lowest pressure occurs downstream close after the front of the cylinder and the highest pressure occurs downstream close after the back of the cylinder. This pressure distribution is due to the fact that the fluid thickness is smaller at the front of the cylinder than at the back, leading to higher flow velocities at the front and, consequently, a lower pressure. perpendicular to the displacement direction. These forces are proportional to the cylinder displacement through the coefficient and , respectively.  Figure 19 shows the evolution of and as a function of , calculated when 0 Hz. In this case, the cylinder is also displaced, according to the scaled by ℎ, and it does not oscillate but only rotates. decreases and increases proportionally to . It can be seen that and evolve similarly to most fluid forces, which usually increase quadratically with the flow velocity, when assuming that the flow velocity is equivalent to .

Added Modal Coefficients under Rotating Conditions
is negative because the associated force tends to move the cylinder away from its equilibrium position.  Figure 19 shows the evolution of and as a function of , calculated when 0 Hz. In this case, the cylinder is also displaced, according to the scaled by ℎ, and it does not oscillate but only rotates. decreases and increases proportionally to . It can be seen that and evolve similarly to most fluid forces, which usually increase quadratically with the flow velocity, when assuming that the flow velocity is equivalent to .
is negative because the associated force tends to move the cylinder away from its equilibrium position. This pressure distribution results in a force towards the minimum fluid thickness, pushing the cylinder away from its equilibrium position, and in another force perpendicular to the displacement direction. These forces are proportional to the cylinder displacement through the coefficient K f and B f , respectively. Figure 19 shows the evolution of K f and B f as a function of Ω, calculated when f = 0 Hz. In this case, the cylinder is also displaced, according to the Φ f1 scaled by h, and it does not oscillate but only rotates. K f decreases and B f increases proportionally to Ω 2 . It can be seen that K f q 1 and B f q 1 evolve similarly to most fluid forces, which usually increase quadratically with the flow velocity, when assuming that the flow velocity is equivalent to Ω. K f is negative because the associated force tends to move the cylinder away from its equilibrium position. Figure 20 shows the normalized values of K f and B f , K f * and B f * , as a function of Ω, calculated when f = 0 Hz based on: where d g is the gap cylinder-wall, and R i is the radii of the cylinder.
where is the gap cylinder-wall, and is the radii of the cylinder. It can be seen that they tend to a constant value when is increased. At lower , * slightly decreases and * increases when is decreased because and are not completely proportional to . At high values of then these coefficients are fully proportional to .   Figure 21 shows the evolutions of , , and as a function of for equal to 0, 0.5, and 1.25 Hz. Similarly, Figure 22 shows their evolutions as a function of when is 0.78, 2.25, and 7 Hz.
As long as is higher than 1 Hz, the cylinder always presents the same value of regardless of and (see Figures 21a and 22a). Consequently, these results suggest that is independent of the state of the flow, as also observed by previous authors [10,[22][23][24], and it also has to be independent of the rotation, based on the present results. may only depend on the distance to the boundaries and the geometry. Exceptionally, at low , increases as presented in Section 3.2 and also observed by [25].
where is the gap cylinder-wall, and is the radii of the cylinder. It can be seen that they tend to a constant value when is increased. At lower , * slightly decreases and * increases when is decreased because and are not completely proportional to . At high values of then these coefficients are fully proportional to .  Figure 21 shows the evolutions of , , and as a function of for equal to 0, 0.5, and 1.25 Hz. Similarly, Figure 22 shows their evolutions as a function of when is 0.78, 2.25, and 7 Hz.
As long as is higher than 1 Hz, the cylinder always presents the same value of regardless of and (see Figures 21a and 22a). Consequently, these results suggest that is independent of the state of the flow, as also observed by previous authors [10,[22][23][24], and it also has to be independent of the rotation, based on the present results. may only depend on the distance to the boundaries and the geometry. Exceptionally, at low , increases as presented in Section 3.2 and also observed by [25]. It can be seen that they tend to a constant value when Ω is increased. At lower Ω, K f * slightly decreases and B f * increases when Ω is decreased because K f and B f are not completely proportional to Ω 2 . At high values of Ω then these coefficients are fully proportional to Ω 2 . Figure 21 shows the evolutions of M f , C f , and A f as a function of f for Ω equal to 0, 0.5, and 1.25 Hz. Similarly, Figure 22 shows their evolutions as a function of Ω when f is 0.78, 2.25, and 7 Hz. presents a linear increase when or increase regardless of the values of or (see Figures 21b and 22b). The rotating flow does not alter the trend of , increasing linearly in general as a function of . However, the rotation alters the exact value of , with higher values of resulting in a higher damping. Assuming that the flow velocity passing over a non-rotating structure is equivalent to the flow induced by a rotating structure, these results agree with previous authors who observed that increased linearly when the flow velocity passing over a vibrating non-rotating structure increased [9,21,23,26]. It can be seen that, in general, increases faster when increases than when increases. Figure 23 presents the turbulence kinetic energy on the midplane of the cylindrical water tank when is 7 Hz and is 0.5 and 1.25 Hz at 0 . Although the flow presents a low degree of turbulence, its increment at higher values of will result in higher VFs and a new pressure distribution that may also contribute to the damping increase. is approximately constant regardless of the value of but it increases linearly when is increased (see Figures 21c and 22c).
is thus independent of and it depends linearly on . Since occurs when the cylinder rotates and increases linearly as a function of , which agrees with the evolution of the Magnus force observed by other researchers [27][28][29], the origin of could be related to the Magnus effect.  As long as f is higher than 1 Hz, the cylinder always presents the same value of M f regardless of f and Ω (see Figures 21a and 22a). Consequently, these results suggest that M f is independent of the state of the flow, as also observed by previous authors [10,[22][23][24], and it also has to be independent of the rotation, based on the present results. M f may only depend on the distance to the boundaries and the geometry. Exceptionally, at low f , M f increases as presented in Section 3.2 and also observed by [25].  C f presents a linear increase when f or Ω increase regardless of the values of Ω or f (see Figures 21b and 22b). The rotating flow does not alter the trend of C f , increasing linearly in general as a function of f . However, the rotation alters the exact value of C f , with higher values of Ω resulting in a higher damping. Assuming that the flow velocity passing over a non-rotating structure is equivalent to the flow induced by a rotating structure, these results agree with previous authors who observed that C f increased linearly when the flow velocity passing over a vibrating non-rotating structure increased [9,21,23,26]. It can be seen that, in general, C f increases faster when f increases than when Ω increases. Figure 23 presents the turbulence kinetic energy on the midplane of the cylindrical water tank when f is 7 Hz and Ω is 0.5 and 1.25 Hz at t = 0T. Although the flow presents a low degree of turbulence, its increment at higher values of Ω will result in higher VFs and a new pressure distribution that may also contribute to the damping increase.
these results agree with previous authors who observed that increased linearly when the flow velocity passing over a vibrating non-rotating structure increased [9,21,23,26]. It can be seen that, in general, increases faster when increases than when increases. Figure 23 presents the turbulence kinetic energy on the midplane of the cylindrical water tank when is 7 Hz and is 0.5 and 1.25 Hz at 0 . Although the flow presents a low degree of turbulence, its increment at higher values of will result in higher VFs and a new pressure distribution that may also contribute to the damping increase. is approximately constant regardless of the value of but it increases linearly when is increased (see Figures 21c and 22c).
is thus independent of and it depends linearly on . Since occurs when the cylinder rotates and increases linearly as a function of , which agrees with the evolution of the Magnus force observed by other researchers [27][28][29], the origin of could be related to the Magnus effect.  A f is approximately constant regardless of the value of f but it increases linearly when Ω is increased (see Figures 21c and 22c). A f is thus independent of f and it depends linearly on Ω. Since A f occurs when the cylinder rotates and A f increases linearly as a function of Ω, which agrees with the evolution of the Magnus force observed by other researchers [27][28][29], the origin of A f could be related to the Magnus effect. Figure 24 shows M f * and C f * corresponding to different Ω as a function of f . These evolutions resemble those of cases where the cylinder is non-rotating, with M f * presenting the same values regardless of Ω and C f * showing the same trend at different Ω but with higher values at higher Ω. Similar to the cases when the cylinder does not rotate, the increase of C f * at low values of f may be partially induced due to the increase in the influence of the VFs.  Figure 24 shows * and * corresponding to different as a function of . These evolutions resemble those of cases where the cylinder is non-rotating, with * presenting the same values regardless of and * showing the same trend at different but with higher values at higher . Similar to the cases when the cylinder does not rotate, the increase of * at low values of may be partially induced due to the increase in the influence of the VFs.     Figure 25 shows the normalized value of A f , A f * , as a function of Ω, based on: A f * presents an approximately constant value regardless of Ω and f , confirming that A f is independent of f and increases linearly with the increase in Ω.  Figure 24 shows * and * corresponding to different as a function of . These evolutions resemble those of cases where the cylinder is non-rotating, with * presenting the same values regardless of and * showing the same trend at different but with higher values at higher . Similar to the cases when the cylinder does not rotate, the increase of * at low values of may be partially induced due to the increase in the influence of the VFs.

Conclusions
This article presents a numerical methodology to calculate the added modal coefficients of a submerged cylinder both when it oscillates without rotation and when it rotates with a whirling motion. It must be noted that a simplified whirling motion was forced in the simulations, presenting slight deviations lower than 1% compared to the case of considering the full whirling motion.
The numerical results of M f , A f , and K f present a close agreement with the corresponding experimental ones. The trend of C f as a function of Ω is captured by the numerical model; however, the numerical predictions present higher values than the corresponding experimental measurements. This deviation is caused by the use of a too-large time step and the additional damping added by the mesh deformation.
The numerical model has permitted the study of particular conditions and allowed conclusions to be reached that could not be reproduced and achieved experimentally, such as: -M f is independent of Ω and f , except at low values of f , when M f increases. C f appears to increase linearly with f and with Ω. A f , K f and B f are independent of f but A f increases linearly with Ω. K f decreases and B f increases with Ω 2 . -When M f and C f are normalized, they tend to a constant value with f , and when A f , K f , and B f are normalized, they also tend to a constant value with Ω. - M f and C f tend to increase for low values of f , probably due to the increasing influence of the VFs.