1. 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 acoustic-structural 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.
2. Numerical Methodology
The following presents a numerical methodology to calculate: (i) the added modal mass,
, and the added modal damping,
, of a standing mode shape when the cylinder does not rotate, and (ii) the added modal mass,
, the added modal damping,
, the added modal stiffness,
, and the two new added modal coefficients proposed in Part 1 [
14],
and
, 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,
and
, which have the same natural frequency and are identical but perpendicular to each other.
Under non-rotating conditions,
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,
, 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)).
Under rotating conditions, the whirling mode shapes can be understood as the combination of
and
, oscillating with the same
and at the same
but at 90° out of phase relative to each other. Thus, having
,
and
expressed as presented in Equations (1) to (3),
,
, and
associated with
must be defined as:
In this case, there are two forces: (i)
in the direction associated with
induced by
,
,
,
, and
as well as their corresponding added modal coefficients, and (ii)
in the direction associated with
induced by
,
,
,
, and
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
and the rotation rather than considering the coupling of both
and
, 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,
in the direction associated with
is only induced by
,
, and
as well as their corresponding added modal coefficients, and
in the direction associated with
is induced by
and
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
is approximately equal to the sine amplitude of
plus the absolute value of the cosine amplitude of
. Using the curve-fitting algorithm, the sine coefficient of
can be calculated and is 1.23 N. Similarly, the sine coefficient of
and the absolute value of the cosine coefficient of
can also be calculated and are 0.7563 and
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:
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.
2.1. Numerical Setup
2.1.1. 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
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:
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:
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.
2.1.2. 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:
where
and
are obtained by solving the following transport equations:
the blending function
is expressed as:
and
as:
where
is the distance to the nearest wall,
is the invariant strain rate, and
is defined as:
The turbulence energy is limited through:
where
is defined as:
All constants are computed via:
and have the following empirical values: (i)
, (ii)
, (iii)
, (iv)
, (v)
, (vi)
, (vii)
, (viii)
, (ix)
.
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
, 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,
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
, a numerical result of
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
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.
2.2. 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.