当前位置:首页 >> 教学计划 >>

Modeling of vortex dynamics in the wake of a marine propeller 相当好的论文


Computers & Fluids 73 (2013) 65–79

Contents lists available at SciVerse ScienceDirect

Computers & Fluids
journal homepage: www.elsevier.com/locate/comp?uid

Modeling of vortex dynamics in the wake of a marine propeller
Roberto Muscari a,?, Andrea Di Mascio b, Roberto Verzicco c,d
a

Istituto Nazionale per Studi ed Esperienze di Architettura Navale (CNR-INSEAN), Via di Vallerano 139, 00128 Roma, Italy Istituto per le applicazioni del calcolo ‘‘Mauro Picone’’ (CNR-IAC), Via dei Taurini 19, 00185 Roma, Italy c Department of Industrial Engineering, Università di Roma ‘‘Tor Vergata’’, Via del Politecnico 1, 00133 Roma, Italy d PoF, University of Twente, Drienerlolaan 5, 7522 NB Enschede, The Netherlands
b

a r t i c l e

i n f o

a b s t r a c t
The ?ow past a rotating marine propeller is analyzed with the aim of establishing limits and capabilities and, hence, the ?eld of applicability of different turbulence modeling approaches for this class of problems. To this purpose the eddy viscosity model of Spalart and Allmaras (1994) [1] and the DES approach [2] have been used. It is shown that the RANSE method can give a very good prediction of global quantities such as thrust and torque, with a relatively small number of grid points. However, when the unsteady ?uctuation of the ?ow or instability processes in the wake are of interest (for noise assessment, for instance), RANSE modeling proves to be too dissipative, as it smoothes out most of the ?nest ?ow features. On the contrary, DES modeling can track the vorticity ?eld for a longer distance and successfully predicts the onset of instabilities in the wake, with excellent agreement with experiments. ? 2012 Elsevier Ltd. All rights reserved.

Article history: Received 9 August 2012 Received in revised form 14 November 2012 Accepted 12 December 2012 Available online 27 December 2012 Keywords: Marine propeller RANSE Detached eddy simulation Vortex instability

1. Introduction The analysis of the ?ow in the wake of a marine propeller is of paramount importance for naval engineering application, it being directly related to vibrations, noise and propulsion performances. Its correct understanding is therefore propaedeutic to the improvement of any existing design or to the development of new concept propellers. The ?ow is very complex [see e.g. 3,4], not only for realistic propeller-hull con?gurations, but also limiting the observation to the simple case of an isolated propeller in a uniform ?ow (open water conditions): in fact, it consists of some principal vortical structures, that can be identi?ed into two vortex systems generated by the root and the tip of the blade (usually very intense) plus a sheet of trailing vorticity produced by the nonconstant circulation along the blade spawn; in addition, a strong vortex emanates from the hub. These vortices interact with each other, giving rise to ?ow instabilities whose proper description is still unclear and under investigation by many research groups. A comprehensive description of the state of the art can be found in Felli et al. [5], where the same ?ow is investigated in a water channel, by ?ow measurements and visualizations. In particular, those authors studied the mechanisms that trigger the instability of the wake by using some reference marine propellers with the same hub and blade pro?les, but different number of blades (2, 3 and
? Corresponding author. Tel.: +39 06 0650299316; fax: +39 065070619.
E-mail addresses: roberto.muscari@cnr.it (R. Muscari), a.dimascio@iac.cnr.it (A. Di Mascio), verzicco@uniroma2.it (R. Verzicco). 0045-7930/$ - see front matter ? 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.comp?uid.2012.12.003

4) and blade loading (advance coef?cient) later de?ned. This investigation allowed to highlight the dependence of the vortex pairing and grouping on the mutual vortex distance. Moreover, they gave an overall description of the way the tip vortices interact with each other and with the hub vortex, when varying the operating conditions. The paper is also very interesting since it relates vortex pairing and grouping to the measured power spectral density, and it shows the way in which the frequency content of the signal is connected to vortex merging and, more in general, to ?ow instabilities. In this paper we investigate on the capability of numerical simulation to reproduce such a complex ?ow. In [6] we considered a ?ow domain much more complex that the one discussed here (a propeller operating in the wake of a fully appended hull) and in that case our experience was that a correctly re?ned RANSE simulation, using the one-equation Spalart and Allmaras [1] eddy-viscosity model, was able to give an excellent description of the ?ow, in terms of both velocity prediction and force coef?cients. However, the physical domain of interest was limited to the zone around the rudder and it was not suf?ciently extended downstream of the propeller plane in order to give a valid estimation on the capability of the RANSE to track the vortex ?laments. Moreover, the propeller operated near its design conditions and the tip vortices showed consistently no sign of instability. A further aspect that triggered the interest for the present analysis is the possibility to perform noise predictions based, for instance, on the FfocwsWilliam Hawking formulation [see 7] where the hydrodynamic

66

R. Muscari et al. / Computers & Fluids 73 (2013) 65–79

data are to be used as input for the noise computations. In this case, we expect the RANSE modeling, by ?ltering all turbulent ?uctuations, to yield spectra with a poor frequency content, which is of course a major problem for noise analysis. Of course, given the ?ow Reynolds number, the use of Direct Numerical Simulation is indisputably out of reach, even on account of the capabilities of new computers in the predictable future. Unfortunately, also a wall-resolved Large Eddy Simulation seems unpractical owing to the discretization requirements in the boundary layer region. Consequently, the only feasible approach to this problem seems to be the Detached Eddy Simulation (DES) [2], where any attempt to directly capture the unsteady dynamics of the ?ow close to rigid walls is replaced by RANSE modeling, while the direct description of the ?ow in the wake (the detached vortices) is retained. The aim of the present paper is to establish limits and capabilities of both the eddy viscosity model by Spalart and Allmaras [1], that proved to be robust and effective for the simulation of the ?ow around a complex geometry (see, for instance, Muscari et al. [6] and of its DES version [2] which should be able to better reproduce the main unsteady features of the turbulent structures in the wake of the propeller. It will be shown that, with a reasonable computational effort, the DES approach can capture the main features of the vortical structures dynamics in the propeller wake, and therefore allows to get the information required, for instance, for the computation of the source terms needed in noise prediction. Moreover, it will be shown that, if considering the running averages of the ?ow ?eld in comparison with the analogous prediction obtained by a RANSE simulation, also the averaged quantities predicted by the RANSE approach are unsatisfactory. On the other hand, when only forces and moments are of interest, the inexpensive (in terms of CPU time) RANSE approach gives practically the same quality of prediction; this suggests that when only forces are required (for instance, for manoeuvrability simulations), a RANSE approach can be conveniently used. 2. Mathematical and numerical models 2.1. Mathematical models The numerical simulation of the ?ow ?eld is performed by the integration of the Reynolds Averaged Navier–Stokes equations written in non-dimensional form in the frame of reference ?xed to the propeller, in terms of absolute velocity components:

In the code, the simulation of turbulent ?ows can be performed by a variety of eddy viscosity models, ranging from the algebraic model by Baldwin and Lomax [8] to the more complex two-equations differential models by Lam and Bremhorst [9] and Chang et al. [10], or by the Detached-Eddy Simulation (DES) or delayed DES (DDES) approaches described in [11]. In all the simulations shown in the paper, the one-equation model by Spalart and Allmaras [1] was used for RANSE simulations, whereas the DES approach [11] was applied for the direct simulation of the larger turbulent structures in the wake. For the sake of completeness, the two models used in the present paper are shortly summarized. The reader is referred to the original papers for a complete description and details. In the one-equation RANSE model [1], the turbulent viscosity mt ~, related to the is computed by means of an auxiliary variable m former by an algebraic equation

~fv 1 with f v 1 ? mt ? m

v3 v ? Cv1
3 3

where

v?

~ m m

?4?

C3 v 1 being a constant. This auxiliary variable evolve according to a convection–diffusion equation with production and destruction source terms

o ~ ~ 1n @m @m ~? ? C b2 jrm ~ j2 ? C b1 ?1 ? f t 2 ? ~?rm ? r ? ??m ? m ? uj @t @ xj r   2 ~ C b1 m ? ft1 DU 2 ? C w1 fw ? 2 ft2 j d

?5?

8 ui <@ ?0 @ xi : @ ui ? ijk xj uk ?
@t @ @ xj @p ??uj ? jlm xl r m ?ui ? ? @ ? xi @ sij @ xj

here, the ?rst term in curly bracket on the right hand side represents diffusion, while the second term in bracket induces propagation of the turbulent front in laminar regions. The other two terms represent production and destruction of turbulence, whereas the last term (called the ‘‘trip’’ term) is used to trigger turbulence at speci?c locations. In all the simulations performed, this last term was always switched off. Finally, the Cs and r are constants, j = 0.41 is the Karman constant, whereas the fs are functions e depending on the distance from the ? wall d and S is another function p??????????? depending on d and on S ? Xij Xij , Xij being the rotation tensor. On the basis of this model, to be used in a RANSE approach, in Spalart et al. [2] a new hybrid model was developed with the aim of achieving an LES simulation of the larger turbulent structures far from rigid boundaries, at the same time overcoming the dif?culties encountered with the use of LES in near-wall regions (Detached Eddy Simulation). This model retains the formal expression of the original model, but the distance function is replaced by

?0

i ? 1; 2; 3 ? 1?

~ ? min?d; C DES D? d

?6?

CDES being another constant and D the larger cell size along the three coordinate curves. 2.2. Numerical models The numerical algorithm is based on a ?nite volume technique with pressure and velocity co-located at cell centers. For the viscous terms a standard second order centered scheme is adopted, whereas for the inviscid part the code has a library of different discretization schemes, a second order ENO (Essentially Non-Oscillatory) scheme, a third order upwind scheme [12], a third order WENO scheme [13] and a fourth order centered scheme [14]. The ENO and WENO algorithms are able to handle ?ows with strong gradients and even discontinuities with remarkable robustness. Nevertheless, they often introduce too much numerical diffusion that can overshadow the physical viscosity, and the number of points required to capture the generation and evolution of some ?ow details can become prohibitive. In these cases, the use of a

where ui are the absolute velocity component in the relative system, xi is the (constant) angular velocity, ri ? xi ? x0 i is the position of point xi with respect to x0 i ; p is the pressure, sij = mT(@ ui/@ xj + @ uj/ @ xi) the stress tensor, mT = 1/Rn + mt, mt being the turbulent viscosity calculated by means of a proper turbulence model (later described) and Rn the Reynolds number (later de?ned). At the physical and the computational boundaries proper conditions are enforced, namely: Solid walls the velocity of the ?uid is set equal to the local wall velocity, i.e.

ui ? ijk xj xk ;

i ? 1; 2; 3

? 2?

In?ow the velocity is set to the value of the undisturbed ?ow

ui ? u1 i ;

i ? 1; 2; 3

? 3?

Out?ow the velocity is extrapolated from inner points.

R. Muscari et al. / Computers & Fluids 73 (2013) 65–79

67

higher order scheme provides the most convenient solution, and, for this reason, all the numerical results shown in this work were obtained with the third order upwind scheme for steady RANSE simulation and the centered fourth-order scheme for the DES simulations. Of course, the main drawback of both schemes is that they lack the non-oscillatory property of the ENO and WENO schemes. However, as no discontinuities in the solution are expected, this is a minor problem; moreover, they seem not to degrade the convergence properties in the sub-iteration for the solution of the implicit problem at each time step (see later on). It should be stressed that the formal spatial accuracy of this method is second-order (owing to the treatment of the viscous terms); nevertheless, the use of a less dissipative scheme for the convective terms allows to follow the formation and convection of strong gradients (in particular, vortex cores) for a considerably longer distance than with the use of a second order-scheme (see, for example [6]). The physical time-derivatives in the governing equations are approximated by a second-order accurate, three-point backward ?nite difference formula [14]; the discrete equation being fully implicit, no stability restriction exists on the time step, that can be chosen on the basis of accuracy alone. In order to obtain a divergence free velocity ?eld at every physical-time step, a dual- or pseudo-time derivative is introduced in the discrete system of equations [15]. The convergence ratio for the inner iteration is accelerated by means of local time stepping, an implicit Euler scheme with approximate factorization [16] and a multi-grid technique [17]. In order to cope with complex geometries or bodies in relative motion, the numerical algorithms were discretized on a block structured grid with partial overlapping, possibly in relative motion. This approach renders domain discretization and grid quality control much easier than with analogous discretization techniques implemented on structured meshes with abutting blocks. Of course, grid connections and overlapping are not trivial, as with standard multi-block approaches, but must be computed in advance. The algorithm used for this purpose is summarized in the following; the interested reader can look for details in [6].  For each internal grid point, we check whether another block is overlapped to the same position. If so, of all possible overlapped grids, we retain only the local ?nest one, and we interpolate the solution on all the other coarser cells, marked as holes.  For each block boundary that has not a simple connection with other blocks, or is not a physical boundary, as before we look for the smallest of all possible donor cells among the other blocks.  Once a cell is identi?ed as a hole or a chimera boundary and the main donor is found, a convex set of eight donors is searched for, and the solution is transferred to the receiver by a tri-linear interpolation. When dealing with unsteady problems with moving boundaries, grid topology has to be recomputed at each time step. In order to render this computation ef?cient, donors and receivers are recomputed only when the blocks they belong to are in relative motion; this means that, for each grid point, we look for possible new donors and receivers only among the cells of blocks that are moving with the one under investigation, at the same time keeping ?xed the connections between any block couple whose relative position does not change in time. Moreover, a nested search algorithm was developed, in which the multigrid structure of the code is fully exploited: 1. We start from the coarsest grid level where we check, for the possible receiver under analysis, every single cell of all other blocks in relative motion as possible donor: this is done in order

to prevent stagnation of the search algorithm described at the next point when applied to non-convex blocks. This search is very quick, because each time that we go down one grid level, the number of cells is reduced by a factor eight (for instance, with three grid levels, the number of grid point to check is reduced by a factor 83 = 512). 2. Once the donor is found on the coarsest grid, this point is used as the starting point of a minimization algorithm along the coordinate lines. Normally, 2 or 3 comparisons of squared distances are enough to locate the closest point in the Cartesian space. This location is then chosen as the initial guess for the next ?ner grid. 3. Once the ?nest level is reached and the closest point is located in the Cartesian space, the search is completed in contravariant coordinates to properly identify the closest point in the computational space and the convex set of eight donors. There is another peculiarity in the algorithm that must be underlined, because it is not standard in chimera-type algorithms. In our approach, the interpolated solution is transferred among different chimera subgrids in a body-force fashion, without the need of fringe points. In particular, for the cells marked as ‘‘hole’’, a source term is added to the discrete equations:

h i j n n ?1 n qn qchimera ? qn chimera ? qchimera ? Dt R ? interp d

?7?

where R is the vector of the residuals, j ? O?10? is a parameter chosen through numerical tests, d is the minimum between (nondimensional) cell diameter and time step, and qn interp is the solution interpolated from the chosen donors belonging to overlapped grids. It can be easily seen that the last term acts as a forcing action that drives the solution toward the interpolated value qn interp . It is very convenient to write the equations in this form because, by doing so, we can retain the data organization of a structured code; consequently, it is very easy to retain, for instance, the multigrid iteration in its standard form, in spite of the holes in the grid. Moreover, these procedures are very convenient when dealing with moving blocks, where some points becomes active after being marked as holes. In fact, for these new points, we must know the solution at previous time steps, in order to compute the time derivatives. If the solution were not forced to the correct value, some special procedure should be masterminded to get the proper value. 3. Test case description The test case chosen for our investigation is the same as that described in Felli et al. [5], where the evolution of the tip and hub vortices in the wake of a marine propeller were experimentally studied in a water tunnel. In that work the authors investigated the effects of blade loading, and therefore of the spiral-to-spiral distance, on the stability of the tip vortices; it was found that the smaller the inter-spiral distance, the sooner instability occurs; moreover, it was found that the evolution of the tip vortices follows a multi-step grouping mechanism. For our purposes, that investigation is particularly convenient: It deals with a relatively simple geometry (a propeller in open water), with a wealth of ?ow visualization details and a variety of different propeller loadings. In fact, since in the numerical simulation the extension of the resolved ?eld past the propeller had necessarily to be limited, we focused our computations on low advance coef?cients (high loading) for a four blade propeller, in order to limit our simulation to those cases where the vortices break-down ‘‘close’’ to the propeller. The propeller geometry is the INSEAN E779A model, i.e. a four blade, ?xed-pitch, right-handed propeller characterized by a nominally constant pitch distribution and a very low skew. Views of the propeller are shown in Fig. 1 whereas the main geometrical

68

R. Muscari et al. / Computers & Fluids 73 (2013) 65–79

Fig. 1. Front and side views of the propeller model.

Table 1 Propeller parameters. INSEAN E779A model Diameter Number of blades Pitch ratio Rake Expanded area ratio Hub ratio D = 0.227 m Z=4 P/D = 1.1 4°350 (forward) 0.689 0.200

features are reported in Table 1. More information on the propeller, can be found in Salvatore et al. [18].

4. Numerical simulation set-up In the numerical simulations, the rotational speed of the propeller has been kept ?xed to a value of n = 25rps and the different advance coef?cients J = U1/nD are obtained by changing the in?ow velocity U1. Unless otherwise speci?ed, all quantities have been cast in non-dimensional form by using as reference values the radius of the propeller (Lref = 0.1135 m), the velocity of the tips of

the blades (Uref = npD ’ 17.85 m/s) and the density of ?uid (qref = 1000 kg/m3). Therefore, the non-dimensional period of revolution is T = 2 p. The Reynolds number is the same of the physical experiment, i.e. Rn = UrefLref/m = 1.78 ? 106. In order to give an estimation of the numerical uncertainty, two grid levels were used, the coarser being obtained by removing every other point from the ?ner. The computation being performed in the rotating non-inertial frame of reference, a steady solution was obtained for the RANSE simulation; for the DES simulation, the solution was computed by including the physical time derivatives, and the chosen time step was always corresponding to a rotation of one degree for the propeller (dt = 2p/360 ’ 0.01745), in the case of the ?ner grid. Fig. 2 shows the building blocks of the mesh: to each blade and to both the fore and the aft parts of the hub, ‘‘O-grid’’ topology was chosen. The rest of the mesh is made up by concentric toroidal blocks that cover the whole computational domain. In order to give an idea of the cell clustering, in Fig. 3 two slices in the planes y = 0 and x = 0 are reported. In particular, the resolved part (?ner grid) of the computational domain extends up to 4.4 radii downstream of the reference plane (x = 0) of the propeller, and this is the region within which we can expect to be able to properly capture the main ?ow features (tip and hub vortices).

Fig. 2. Grid topology: overview (left) and zoom of the near ?eld (right).

R. Muscari et al. / Computers & Fluids 73 (2013) 65–79

69

Fig. 3. Details of the volume mesh. Section y = 0 (left) and x = 0 (right).

Table 2 Distribution of computational cells among the grid blocks. Four propeller blades Propeller hub Near ?eld (tip vortices) Near ?eld (buffer zone) Near ?eld (all the rest) Far ?eld (background) Total 2.36 M 0.64 M 3.11 M 2.13 M 1.89 M 1.18 M 11.31 M (20.9%) (5.6%) (27.5%) (18.8%) (16.7%) (10.4%)

Therefore, we evaluated the grid uncertainty presented in the next sections, following Roache [19]. In particular, we ?rst estimated the solution variation as

E?

f2 ? f1 1 ? rp

?8?

About 32 cells were put in the boundary layer, the ?rst point being set at a distance from wall such that y+ < 1 in wall unit. Apart from the blocks ?tted around the geometry, most of the cells are placed in the toroidal block that is used to track the tip vortices and in the following; another external toroidal block is used as a buffer between the ?ner mesh near the propeller and the outer, coarser, background mesh. The different blocks sum up to a total of 11 M cells and their relative weight is detailed in Table 2. Fig. 4 shows the mesh density in a zone encompassing the ?rst two vortex cores downstream the propeller for J = 0.38 and J = 0.71 on the plane y = 0. In the same ?gure, the iso-levels report the pressure distribution and the vectors correspond to (u ? U1, 0, w). Despite the dif?culty to de?ne the extension of a vortex core rigorously, it can be seen that the width of each ?lament is described by about 10–12 cells. Finally, for the initial conditions, the steady RANSE simulations started from a uniform ?ow condition (with U1 prescribed by the advance coef?cient under analysis), whereas the unsteady DES simulations started from the corresponding (i.e. with the same advance coef?cient) RANSE computation. As already noted, the RANSE simulation is a steady one and consequently the computational cost is relatively small: on a cluster of 128 AMD Opteron processors it took about 20 ‘‘wall clock’’ hours in order to reach iterative convergence on the ?nest grid. Conversely, the DES were unsteady simulations and were carried on for about 12 propeller revolutions. On the same cluster as before a typical DES lasted about 350 ‘‘wall clock’’ hours.

where f1 and f2 are, respectively, the ?ne-grid and coarse-grid solutions, r the re?nement factor (r = 2, as stated before) and p is the formal order of accuracy (p = 2, for the present code). Then, we evaluated the uncertainty by

U N ? F s jEj

?9?

where Fs is a safety factor that we take, according to Roache, equal to Fs = 3. As to iterative convergence for each sub-step, it was found to be always negligible with respect to grid uncertainty. 5. Results 5.1. Thrust and torque As a ?rst check of the numerical results and in order to compare different approaches for the prediction of the ?ow around the propeller, we considered the open water characteristics of the model. In Fig. 5 several numerical results for different values of the advance coef?cient are superposed to the experimental curves [18]. Concerning the thrust, the predictions obtained with the different approaches (RANSE on both coarse and ?ne meshes, DES on the ?ne mesh) are essentially identical for all values of the advance coef?cient. It can be observed from the ?gure that the difference between computed and measured values is almost constant in the whole simulation range, although it is smaller in percent of the computed values at low values of J (being about 4% of the experimental value) and grows with J up to a difference of about 15% for J = 0.88, where both the computed and the measured values of the force become very small. In addition, it has to be noted the lack of uncertainty estimation of the experimental data. The torque shows a greater sensitivity to mesh re?nement, its value being strongly dependent on a correct estimation of the viscous forces in the boundary layer. In particular, it can be noticed a shift of the whole curve towards lower values as the grid gets ?ner. In any case, the experimental data fall within the range of numerical uncertainty. The use of a DES model, in place of the RANSE simulation, produces some effects only for lower values of J, the differences being irrelevant for higher values. The agreement among the computed forces for the different turbulence modeling (DES and RANSE) is not surprising, since in both cases the boundary layers have been modeled using Spalart

4.1. Uncertainty evaluation In order to get a rigorous evaluation of the numerical uncertainty at least three mesh levels would be necessary, but for the present work this was unfeasible. In fact, a further reduction of the coarse mesh would have lead to a too coarsened grid, unable to resolve many important ?ow features, whereas further re?ning of the ?ne mesh would have required too large computing resources. p??? On the other hand, grid coarsening with noninteger factors (say 2) would yield blocks with a number of points un?t for multigrid algorithm (it must be a number m ? 2n, with m and n both integers).

70

R. Muscari et al. / Computers & Fluids 73 (2013) 65–79

Fig. 4. Mesh density in the core of the tip vortices. J = 0.71 (up) and J = 0.38 (down).

Fig. 6. Vortical structures visualization for J = 0.71.

Fig. 5. Measured and predicted open water characteristics of the E779A propeller.

and Allmaras [1], Spalart et al. [2], and therefore the observed small differences are to be ascribed only to the different vorticity ?eld in the wake of the propeller. It can be concluded that, from the point of view of global quantities estimation, there seems be no point in using a more accurate turbulence model or even a highly re?ned mesh. It will be shown that the situation is strongly different when considering the wake structure, turbulence dynamics and pressure ?uctuations. 5.2. Flow ?eld in the wake The analysis of the ?ow ?eld has been carried out for four values of the advance coef?cient, namely J = 0.71, 0.45, 0.38, 0.2. The ?rst

two values have been considered in order to try to reproduce some of the ?ndings of Felli et al. [5] and to validate the numerical simulations. In particular, for the four-blade con?guration at J = 0.71, the detailed frequency analysis performed in the experimental work illustrates the process of energy transfer from the blade harmonic to the shaft harmonic along the longitudinal direction, owing to vortex grouping. Unfortunately this process takes place after a long distance downstream of the propeller (in the original reference it develops at a distance x = 12.65R), whereas the re?ned computational mesh terminates only at x = 4.4R. In order to be able to capture at least the ?rst stages of the vortex coupling phenomenon described in the experimental work, we repeated the simulation for J = 0.45 (which is the lowest value of the advance coef?cient reported in Felli et al. [5], and for J = 0.38. The last considered value J = 0.20 has been simulated, despite the lack of measurements for validation, as it better illustrates the need to perform a DES

R. Muscari et al. / Computers & Fluids 73 (2013) 65–79

71

Fig. 7. Turbulent kinetic energy for J = 0.71; resolved (left), modeled (according to Eq. (10)) (right).

simulation instead of more classical RANSE calculation for high blade loading (off-design) conditions. 5.2.1. J = 0.71 A general overview of the vorticity ?eld in the wake of the propeller is given in Fig. 6. The surface shown in the ?gure corresponds to the non-dimensional value k2 = ?2 of the second largest invariant of S2 + X2 (S and X being the symmetric and antisymmetric components of ru, respectively [20]) and it is colored1 by pressure levels. After a slight, initial reduction of the radial position of the vortex cores, caused by the acceleration of the ?ow behind the propeller, the helices formed by the tip vortices remain located on a cylinder of almost constant radius up to the end of the re?nement blocks. This simulation reproduces the experimental ?ndings visualized in Felli et al. [5] for J = 0.65 and J = 0.75, where the destabilization process takes place beyond the limit of ?nely discretized domain in the numerical simulation. In order to have at least a qualitative idea of the level of resolution obtained with the grid employed in the simulations, we computed the amount of resolved kinetic energy (that is, TKEres ? 1 ?u02 ? v 02 ? w02 ?) in the wake of the propeller and com2 pared it with a quantity that can give an estimation of the modeled part (as well known, in a Boussinesq’s approximation of turbulent stresses, the modeled kinetic energy is absorbed in the pressure terms and cannot be discerned from it)

Fig. 8. Probes positions relative to the PSD in Fig. 9.

TKEmod ?

1 2

@ u @ v @ w mturb @x ? @y ? @z :













?10?

These two quantities are compared in Fig. 7, where it can be seen that the modeled part of the TKE is reasonably small when compared to TKEres. This shows that cell size is small enough to resolve the largest velocity scale ?uctuations. It has to be noted that the TKEres has also a low maximum value (with respect to the cases reported in the next sections) which corresponds to a very stable wake. A partial validation of the numerical results can be done by comparing the computed power spectral density of the kinetic energy with the analogous experimental data reported in Felli et al. [5]. Numerical data were collected as in the experiments, i.e. we set a series of ‘‘probes’’ in the inertial frame of reference at a distance r = 0.9R from the axis of the propeller and at several stations downstream (Fig. 8) and the time histories of the kinetic energy were recorded; then, we have computed the power spectral density (PSD) of the kinetic energy K = 0.5(u2 + v2 + w2). The result for three probes are shown in Fig. 9, where the wave numbers have been scaled so that the blade frequency is j = 1. In particular, the
1 For interpretation of color in Figs. 6, 7, 9, 16 and 22, the reader is referred to the web version of this article.

?rst two probes have been selected in order to make a direct comparison with experiments (see Fig. 40 in Felli et al. [5] and the third one is placed at the farther distance allowed by the numerical grid, being near the boundary of the inner ?ne computational mesh. Several observations can be done from the reported graphs. For X = 0.35R, very close to the blade tips, the spectrum exhibits sharp peaks at the blade frequency and all its multiples, indicating a low level of noise superimposed to the tonal components. Although the experimental spectra in Felli et al. [5] are truncated at a frequency only slightly higher than twice the blade frequency, at this probe position they also show well de?ned peaks for j = 1 and j = 2. If we inspect the data at X = 3.55R, the energy transfer thoroughly described in Felli et al. [5] is beginning to appear also in the numerical data; in particular, a smaller peak for j = 0.5 is evident in the spectrum indicating an incipient pairing between the tip vortices. The main difference between CFD and measurements is the presence (in the experiments) of a strong peak for j = 1.5 which is visible in the numerical spectrum although not as strong as in the experiments. This is probably to be ascribed to the perfectly uniform in?ow in the numerical simulations, that delays the instability onset; in fact, as shown for the case J = 0.45 in the following, even a small level of turbulence can trigger vortex coupling and instability closer to the propeller. The last station that could be numerically considered, X = 4.0R, has no counterpart in Felli et al. [5], but it is very interesting as it shows both a more pronounced peak at j = 1.5 and a ?rst appearance of a peak at the shaft frequency (j = 0.25) which is, according to the experiment, the only contribution to the spectrum when the process of energy transfer completes. In Fig. 9 is also shown, through the red line j?0.9, that the power-law decay found in the experiments is very well captured. As to the question regarding the need of a (very CPU expensive) DES simulation with respect to a much faster RANSE approach, a comparison between the predictions of the ?ow ?elds, obtained by the two different turbulence models is given in Fig. 10; here, the results from the steady RANSE computation are juxtaposed to the ?eld obtained by averaging 3600 instantaneous ?elds (i.e. 10

72

R. Muscari et al. / Computers & Fluids 73 (2013) 65–79

Fig. 9. J = 0.71 – Power spectral density of the kinetic energy at different stations downstream the propeller.

Fig. 10. J = 0.71. Top: section x–z of the three-dimensional axial velocity ?eld; bottom: vorticity ?eld visualization. RANSE (left) and running averages from DES (right).

revolutions) from the unsteady DES. In particular, the top frames show a longitudinal cut of the axial velocity ?eld. The main ?ow characteristics are the same for the two approaches, with a well de?ned zone of accelerated ?ow, in which the traces of the tip vortices and the wakes of the blades, together with a central zone past the hub of slower ?uid, are evident. However, all these features are very smoothed in the RANSE simulation, and some are completely missing (like the very low velocity zone in the core of the hub vortex). Moreover, when comparing the resolved vortical structures (obtained by the same method as in Fig. 6), it is clear that the RANSE model yields tip vortices than vanish after a length of about a propeller diameter. In our opinion, these differences are caused by a too high level of viscosity introduced by the RANSE turbulence model; in order to verify this hypothesis, we have plotted in Fig. 11 the total viscosity (that is, mT = 1/Rn + mt) and the y-component of the vorticity in the plane x–z for both the DES and the RANSE. The RANSE simulation predicts a small mT near the tip of the blade, but it rapidly grows reaching the maximum value about one spiral downstream. Also, the area around the vortex core where the mT has a substantial value widens along the ?lament. On the contrary,

the DES produces a higher total viscosity near the tip but already after a quarter of spiral the level of viscosity has reduced and the area affected by it has reached its de?nitive extension. Consistently, the level of vorticity, which is very similar near the propeller for the two simulations, drops very fast for the RANSE whereas has a slower decay for the DES.

5.2.2. J = 0.45 This value of the advance coef?cient is the lowest one for which visualizations are reported in Felli et al. [5], and has been analyzed in order to check whether the numerical simulation was able to capture the vortex pairing phenomenon, which for this J begins (by a mere observations of the photographs) around X = 0.3R. In Fig. 12, similarly to the case J = 0.71, the comparison between RANSE and DES simulation is reported. In particular, the bottomright picture is the analogous of the upper-right frame of Fig. 8 in Felli et al. [5], where the pairing of tip vortices and the incipient instability of the hub vortex are visualized by inducing cavitation of the vortex cores (by lowering the ambient pressure).

R. Muscari et al. / Computers & Fluids 73 (2013) 65–79

73

Fig. 11. J = 0.71 – section x–z. Total viscosity (upper half) and y-component of vorticity (lower half). Top: DES; bottom: RANSE.

The higher thrust obtained for this value of the advance coef?cient J (and therefore the higher blade loading) manifests both in the faster deformation of the wakes that, in the inner part of the slipstream, become rapidly aligned to the ?ow, and in the stronger tip vortices (it should be noted that the velocity contour levels are the same as in Fig. 10). In comparison with the DES, the RANSE approach seems to introduce even more viscosity than in the case

J = 0.71, as in the last part of the resolved region of the wake the velocity gradients, still very strong in the DES, are completely blurred. It is also signi?cant that, despite the fact that the tip vortices are stronger than in the previous case, they can be tracked for a shorter distance (see the bottom-left frame of Fig. 10 for comparison) before being dissipated. This trend is con?rmed by the following cases, J = 0.38 and J = 0.2, and evidently means that the

Fig. 12. J = 0.45. Top: section x–z of axial velocity; bottom: vorticity ?eld visualization. RANSE (left) and running averages from DES (right).

74

R. Muscari et al. / Computers & Fluids 73 (2013) 65–79

Fig. 13. J = 0.45. Vorticity ?eld obtained by adding a random velocity ?eld to the cells at x < ?1.

growth of turbulent viscosity in the Spalart and Allmaras [1] model triggered by the velocity gradients is, for this kind of computation, de?nitely too high. On the other hand, the DES permits to track the vortices up to the end of the resolved mesh, and it effectively shows both the onset of the vortices instability and the start of the pairing process described in the experiments. A visual comparison of the bottom-right part of Fig. 12 with the corresponding photograph in Felli et al. [5] shows a striking similarity, with the features of the con?guration marked by red ellipses identical in both data sets. Actually, there is a difference in the starting point of the transition to instability, the numerical one being retarded with respect to the experimental one. In order to explain this seemingly stronger stability of the ?ow in the numerical calculation, even with DES modeling, we made some additional simulation to check whether the instability onset could be affected by disturbances in the incoming ?ow. In fact, the above simulations report an ideal situation with a perfect uniform in?ow, whereas the experiments are conducted in a water channel where the incoming ?ow has necessarily both a background turbulence level and a certain level of nonuniformity, not to mention the fact that the visualizations were obtained by inducing phase transition (cavitation) in the water. To check this hypothesis, we repeated the previous calculation by adding to the uniform in?ow a random velocity ?eld:

turbulence would probably lead to a better similarity with the experimental results but it goes beyond the scope of this paper. With respect to the case J = 0.71, the intensity of the resolved TKE (Fig. 14) has grown, in particular towards the end of the resolved mesh where the instability process begins, whereas the modeled TKE is limited to lower levels. In the plot of the resolved TKE, it is interesting to note the start of the pairing process between the two couples of tip vortices that is revealed with extreme clarity. This pairing process can also be seen from the analysis of the power spectral density of the kinetic energy. At x = 0.35R (Fig. 15) well de?ned peaks at the blade frequency and its multiples can be observed, as expected because of the small distance from the propeller plane. However, at a distance x = 2.5R most of the peaks have vanished and the whole curve has shifted to higher values, suggesting a higher noise content. The only peaks that still characterize the spectrum are for j = 1 plus, though with a smaller magnitude, j = 0.5 and j = 1.5, these latter indicating the process of incipient vortex coupling already observed previously. At x = 3.55R the whole spectrum has further increased in magnitude and, apart from a small peak corresponding to the blade frequency, no other harmonic component seems to characterize it. As for the case J = 0.71, we found a power-law decay j?0.9. 5.2.3. J = 0.38 For this value of the advance coef?cient no experimental data are available, neither in form of tip vortices visualization nor in form of PSD of the velocity ?eld. However, we ?nd it interesting in that it permits, more clearly than the other cases, to assess the capacity of the numerical method to follow the onset of the instability and the process of vortex coupling in the propeller wake. In order to better distinguish between the two couples of interlacing vortices, in Fig. 16b each pair has been represented by a different color. In this way it is possible to observe the initial stages of the instability mode described in Felli et al. [5] as ‘‘leapfrogging’’ phenomenon (see in particular Fig. 12 of the cited paper, referring to the case J = 0.65): two adjacent vortex ?laments which have initially the same radius (ellipse tagged with ‘‘1’’ in Fig. 16b) in?uence each other so that the aft ?lament is pulled forward, pass inside the other vortex spiral and ?nd itself in the fore position (ellipse tagged with ‘‘2’’). In the experiments the process repeats itself but we could numerically capture only the ?rst occurrence of the phenomenon. It is worth noting that after the instability process starts, the exact position of the vortex ?laments oscillates in time and that the con?guration shown in Fig. 16b is just a snapshot at a precise instant. As a matter of fact, using the time-averaged values of the velocity ?eld for the computation of k2 would lead to Fig. 16a. As for the previous cases, we compared the amount of the resolved kinetic energy to the estimated modeled portion of it

un inflow ? u1 ? 0:2R?t ; x; y; z?jx?x1 ju1 j
where R is a random vector function such that jRj < 1. By this perturbation we have a maximum associated kinetic energy equal to 4% of the undisturbed ?ow, not far from the background value declared for the experimental facility. Of course, this is a very crude approximation of a turbulent ?ow, but it is suf?cient for our purposes. The outcome of the computation is reported in Fig. 13 where it is evident that this small disturbance produces dramatic effects on the instability onset. A better characterization of the incoming

Fig. 14. Turbulent kinetic energy for J = 0.45; resolved (left), modeled (according to Eq. (10)) (right).

R. Muscari et al. / Computers & Fluids 73 (2013) 65–79

75

Fig. 15. J = 0.45 – Power spectral density of the kinetic energy at different stations downstream the propeller.

Fig. 16. Tip vortex visualization for J = 0.38: (a) time-averaged ?eld from DES, (b) snapshot of instantaneous ?eld.

Fig. 17. Turbulent kinetic energy for J = 0.38; resolved (left), modeled (according to Eq. (10)) (right).

(Fig. 17). For this blade loading, the amount of resolved kinetic energy is signi?cantly larger than the modeled part, and its maxima relevantly increase in comparison with the previous test cases. The observation of the resolved kinetic energy ?eld clearly reveals the vortex grouping phenomenon when moving downstream. The power spectral density is reported in Fig. 18. It is clear from the graph that the periodic behavior of the function dominates close to the propeller disk, where the vortices still retain their helicoidal shape. Position X = 1.9R, instead, is the place where vortex coupling begins and the contribution from the component whose frequency is one half of the blade frequency appear. Moreover, at this position still exists a portion of the PSD graph whose decay is very close to the slope that was observed, numerically and experimentally, for the two previous cases.

When moving further downstream (X = 2.5R), the periodic character of the ?ow is completely lost, as dominant components are no longer identi?able. Here the decay of PSD with frequency is much faster than j?0.9. 5.2.4. J = 0.20 The last test case we considered is a condition for which the advance coef?cient is very low, and therefore the blade loading very high. Tip vortex evolution can be observed from the right column of Fig. 19, where the instantaneous ?ow ?eld represented in terms of k is reported together with a longitudinal section of the averaged ?ow ?eld in terms of axial velocity. The results reported on the right of Fig. 19 were computed with the DES modeling and can be compared with analogous simulation performed with the RANSE simulation of the

76

R. Muscari et al. / Computers & Fluids 73 (2013) 65–79

Fig. 18. J = 0.38 – Power spectral density of the kinetic energy at different stations downstream the propeller.

Fig. 19. J = 0.2. Top: section x–z of axial velocity; bottom: instantaneous vorticity ?eld visualization. RANSE (left) and DES (right).

left of the same ?gure. Even more clearly than in the previous cases, these results show that a RANSE simulation, although can quickly yield good predictions of forces and moments acting on the propeller, completely fails to give correct information about the wake evo-

lution, whose details are of paramount importance if the noise generated by the propeller is to be estimated. The comparison of the resolved kinetic energy and the estimation of the modeled part is reported in Fig. 20. For this test case,

R. Muscari et al. / Computers & Fluids 73 (2013) 65–79

77

Fig. 20. Turbulent kinetic energy for J = 0.2; resolved (left), modeled (according to Eq. (10)) (right).

Fig. 21. J = 0.2 – Power spectral density of the kinetic energy at different stations.

the resolved kinetic energy clearly grows and dominates over the modeled part. As to the PSD, the graphs in Fig. 21 show that the periodic character of the signals can be distinguished only for X = 0.35R; at the two following probes X = 0.50R and X = 0.65R, the periodicity of the phenomenon leaves its mark only for the persistence of the component with the blade frequency. Moving further downstream, any trace of forced periodicity almost disappears, and the components with frequency lower than the blade frequency (and in particular the shaft frequency 1/4) begin to grow. This last behavior of the ?ow ?elds corresponds to what is observed in the far ?eld for higher values of the advance coef?cient in the experimental obser-

vation, and its is due to the complete grouping of the four tip vortices and their intertwining with the hub vortex. These phenomena can be highlighted if the PSD are measured in the rotating reference system, because this change of reference completely purges the signals from the forced periodicity due to rotation. In Fig. 23, the spectra are measured at several positions in typical regions of the domain in the rotating system, and then averaged. The regions considered are (see Fig. 22):  a cylindrical zone behind the hub (red zone in the ?gure);  a toroidal region around R = 0.5 (mid-blade) in the wake of the blade (green zone);

78

R. Muscari et al. / Computers & Fluids 73 (2013) 65–79

the nonlinear interaction of two signals, one with the blade frequency j = 1 and another one with the shaft frequency j = 0.25. This shape of the PSD can be more clearly observed in the other two regions, tip-far and tip-near. 6. Conclusions The capabilities of numerical simulations with different turbulence models (RANSE and DES) to predict the complex ?ow past an isolated propeller have been assessed. In order to investigate the wake structure and turbulence dynamics, different operational conditions have been considered, starting from a moderate blade loading, at J = 0.71, up to very high loading at J = 0.2. A direct comparison between the Spalart–Allmaras eddy-viscosity model and the DES approach (resting itself on the Spalart–Allmaras model) shows that the RANSE approach dissipates the vorticity ?eld very quickly and, being all other parameters (the computational mesh and integration scheme) identical to the corresponding DES, this over-dissipation is to be completely ascribed to the eddy-viscosity modeling. Furthermore, the growth of turbulent viscosity, strictly connected to the velocity gradients, is such that the stronger the tip vortices (for low values of J), the sooner they are dissipated (and, hence, the shorter distance they can be tracked). On the contrary, the DES method allows to capture the tip vortices evolution as long as the mesh is reasonably re?ned with a good qualitative and quantitative agreement with experiments. In particular, the initial stages of the instability pattern, with two consecutive vortex ?laments swapping their relative position because of the mutual interaction, can also be reproduced (see cases J = 0.45, 0.38) and follow with striking similarity the ?ow visualizations from water channel experiments. In this respect, the main

Fig. 22. Probe locations for the power spectral densities in the rotating reference system.

 a toroidal region behind the tip of the blade for 0.16R < x < 0.75R (tip-near, blue zone);  a toroidal region behind the tip of the blade for 0.76R < x < 1.28R (tip-far, cyan zone). These probes are obviously ?xed with respect to the blade. It can be seen that in the hub region there is no dominant component and besides two portions for the spectrum with slope j?0.9 and j?5/3 can be identi?ed. In the mid-blade region, instead, components with frequencies around j = 0.75 and j = 1.25 appear, together with their multiples, which is to be probably attributed to

Fig. 23. J = 0.2 – Power spectral density of the kinetic energy at different stations in the rotating frame of reference.

R. Muscari et al. / Computers & Fluids 73 (2013) 65–79

79

limit of the computations seems to be the length of the resolved mesh, that is the number of computational volumes. Through a random disturbance to the velocity ?eld upstream the propeller, it has been shown that the slight delay in the onset of the instability for the numerical simulation of case J = 0.45 can be explained by the fact that the numerical in?ow is perfectly uniform, whereas the in?ow to the propeller in the water channel is obviously affected by some level of turbulence. The analysis of the power spectral densities of the signals at different downstream locations has also been performed. Although the resolved ?eld was not suf?ciently extended and we could provide information only up to X = 4.0R, numerical simulations were able to con?rm the process of energy transfer, from the blade- to the shaft-frequency components, associated to the vortex pairing observed in the ?ow visualizations. Moreover, the same power-law decay measured in the experiments was found in the numerical spectra, as a veri?cation of the correct representation of the energy transfer between turbulent scales. A further check on the quality of the DES is given in terms of comparison between the resolved part and (an estimate of) the modeled part of the turbulent kinetic energy: the latter proves to be smaller than the former when a fully turbulent ?ow is to be expected, con?rming that the grid resolution was enough to resolve the largest scales of velocity ?uctuations. Finally, the results obtained for the case J = 0.2, featuring a vorticity ?eld extremely rich in terms of structures of different sizes, have highlighted that, whenever the frequency content of the signals spectra is important (for instance when the noise generation is to be estimated), the DES approach is the better trade-off between an unfeasible LES/DNS and the too high turbulent viscosity introduced by RANSE. On the other hand, as far as global quantities such as thrust and torque are concerned, the use of DES instead of less demanding RANS is pointless, as both approaches yield very good agreement with experimental data. Acknowledgments This work was partially supported by the Italian Ministry of Education, University and Research through the research project RITMARE. Numerical computations presented here have been performed on the parallel machines of CASPUR Supercomputing Center (Rome); their support is gratefully acknowledged.

References
[1] Spalart PR, Allmaras SR. A one-equation turbulence model for aerodynamic ?ows. La Recherche Aérospatiale 1994;1:5–21. [2] Spalart PR, Jou W-H, Strelets M, Allmaras SR. Comments on the feasibility of LES for wings, and on a hybrid RANS/LES approach. Columbus (OH): Greyden Press; 1997. [3] Kerwin JE. Marine propellers. Annu Rev Fluid Mech 1986;18:367–403. [4] Conlisk AT. Modern helicopter aerodynamics. Annu Rev Fluid Mech 1997;29:515–67. [5] Felli M, Camussi R, Di Felice F. Mechanisms of evolution of the propeller wake in the transition and far ?elds. J Fluid Mech 2011;682:5–53. [6] Muscari R, Felli M, Di Mascio A. Analysis of the ?ow past a fully appended hull with propellers by computational and experimental ?uid dynamics. ASME J Fluids Eng 133 (6). http://dx.doi.org/10.1115/1.4004215. [7] Ianniello S. New perspectives in the use of the Ffowcs Williams–Hawkings equation for aeroacoustic analysis of rotating blades. J Fluid Mech 2007;570:79–127. [8] Baldwin B, Lomax H. Thin-layer approximation and algebraic model for separated turbulent ?ows. In: AIAA, 16th aerospace sciences meeting, Huntsville, Alabama, USA, AIAA Paper 1978-257, 1978. [9] Lam CKG, Bremhorst K. A modi?ed form of the k-epsilon model for predicting wall turbulence. ASME J Fluids Eng 1981;103(3):456–60. [10] Chang KC, Hsieh WD, Chen CS. A modi?ed low-Reynolds-number turbulence model applicable to recirculating ?ow in pipe expansion. ASME J Fluids Eng 1995;117:417–23. [11] Spalart PR. Detached-eddy simulation. Annu Rev Fluid Mech 2009;41:181–202. [12] Van Leer B. Towards the ultimate conservative difference scheme V. A secondorder sequel to Godunov’s method. J Comput Phys 1979;32(1):101–36. [13] Jiang G-S, Shu C-W. Ef?cient implementation of weighted ENO schemes. J Comput Phys 1996;126(1):202–28. [14] Di Mascio A, Broglia R, Muscari R. Prediction of hydrodynamic coef?cients of ship hulls by high-order Godunov-type methods. J Mar Sci Technol 2009;14:19–29. [15] Merkle CL, Athavale M. Time-accurate unsteady incompressible ?ow algorithm based on arti?cially compressibility. AIAA paper 87-1137. [16] Beam RM, Warming RF. An implicit factored scheme for the compressible Navier–Stokes equations. AIAA J 1978;16:393–402. [17] Favini B, Broglia R, Di Mascio A. Multi-grid acceleration of second order ENO schemes from low subsonic to high supersonic ?ows. Int J Numer Method Fluids 1996;23:589–606. [18] Salvatore F, Pereira F, Felli M, Calcagni D, Di Felice F. Description of the INSEAN E779A propeller experimental dataset. INSEAN Tech. rep. 2006-085; 2006. [19] Roache PJ. Quanti?cation of uncertainty in computational ?uid dynamics. Annu Rev Fluid Mech 1997;29:123–60. [20] Jeong J, Hussain F. On the identi?cation of a vortex. J Fluid Mech 1995;285:69–94.


赞助商链接
相关文章:
更多相关标签: