Development of hybrid URANS–LES methods for flow simulation in the ship stern area

Ocean Engineering 38 (2011) 1831–1838

Contents lists available at SciVerse ScienceDirect

Ocean Engineering journal homepage: www.elsevier.com/locate/oceaneng

Development of hybrid URANS–LES methods for flow simulation in the ship stern area N. Kornev a, A. Taranov a, E. Shchukin a,n, L. Kleinsorge b a b

Chair of Modeling and Simulation, University of Rostock, 18057 Rostock, Germany Chair of Shipbuilding, University of Rostock, 18057 Rostock, Germany

a r t i c l e i n f o

abstract

Article history: Received 26 October 2010 Accepted 10 September 2011 Editor-in-Chief: A.I. Incecik Available online 28 September 2011

The paper presents the results of the application of a new hybrid URANS–LES method for the investigations of the ship wake behind the tanker KVLCC2. The switching between URANS and LES models is based on the ratio between the turbulence scale and the cell size of the mesh. Ship resistance, fields of the axial velocity and turbulent kinetic energy in the propeller plane are calculated and compared with measurements. Much attention is paid to the analysis of the unsteady velocities, their PDF distributions and spectra. Numerical analysis shows that the instantaneous velocities deviate substantially from their mean values which are usually used as the estimated velocities in modern engineering methodologies. The thrust variation in the unsteady wake is more than twice as large as that in the time averaged (frozen) wake. The results of the present study point out that the unsteadiness in the wake behind full ships can be very large and should be taken into account when propulsion and unsteady loadings are determined. & 2011 Elsevier Ltd. All rights reserved.

Keywords: Ship wake Propeller URANS LES Hybrid methods

1. Introduction Determination of the velocity field in the ship wake is one of the most important problems of ship hydromechanics since the wake has strong influence on the ship performance. From one side, the wake is responsible for the flow velocity decrease in the propeller disk, which results in the thrust increase and in the improvement of the overall propulsion efficiency. From the other side, a nonuniform velocity field causes variation of the propeller thrust in time and thus strong vibrations in the stern area. Propeller cavitation is another critical phenomenon in ship hydromechanics which is strongly influenced by the wake. Numerical simulation of the wake has attracted the attention of CFD experts for a long time. A substantial success has been achieved in this field in the last two decades. Today, the averaged velocity field can be computed with high accuracy and the discrepancy between the numerics and the measurements is comparable with the tolerance of the experimental equipment. To get overview of the state of the art in this area, the reader is referred to contributions presented in the Gothenburg (2010) and NuTTS (2010) conference.

n

Corresponding author. E-mail addresses: [email protected] (N. Kornev), [email protected] (A. Taranov), [email protected] (E. Shchukin), [email protected] (L. Kleinsorge). 0029-8018/$ - see front matter & 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.oceaneng.2011.09.024

An important feature of the works done so far is the application of URANS (unsteady Reynolds averaged Navier Stokes) method which is capable of capturing steady effects and large scale unsteadiness. However, this technique is not able to resolve small scale flow oscillations due to large diffusivity which is an unavoidable feature of URANS closure models. These oscillations are caused by complicated vortex structures arisen due to flow separations on the hull and shedding of the boundary layer in the stern area. The time averaged vortex structure is well reproduced in URANS calculations, which is confirmed by a good agreement between the numerical simulations and the measurements of the mean velocities. Also, the time averaged fluctuation parameters like Reynolds stresses and the turbulent kinetic energy are predicted relatively well using advanced turbulent models such as the Reynolds stress models, see Kim and Rhee (2002). This allows one to accurately compute the variances of fluctuations, but not their amplitudes or spatial and temporal distributions. It can be a critical point for propeller design since flow, cavitation and thrust depend on the instantaneous spatial distributions of velocity in the incident flow. In present practical design methods the unsteady effects are partly taken into account. The velocity vector u at any point x on the propeller blade depends on the angular position of the blade y which itself depends on time, i.e. u ¼ uðx, yðtÞÞ. Unsteadiness of the propeller flow is purely due to the rotation of the propeller blade through the nonuniform time independent (‘‘frozen’’) wake (see, for instance, Kornev, 2010, Chapter 12). We usually can assume

1832

N. Kornev et al. / Ocean Engineering 38 (2011) 1831–1838

that the incident flow does not explicitly depend on the time. This assumption is quite acceptable for ships with moderate block coefficients since the physical phenomena causing flow oscillations do not play a big role for such ships. However, the case of full-bottomed ships should be reconsidered. Because of the disadvantages of URANS modeling mentioned above, the most promising approach to work with such ships is large eddy simulation (LES), which is already widely used for research purposes. On the contrary to other engineering fields, typical Reynolds numbers in ship hydromechanics are very large even at model scales. The grid resolution necessary for a pure LES is so huge that it makes the direct application of LES impossible. A practical solution of this problem is the use of hybrid URANS–LES methods, where the near body flow region is treated using URANS and far flow regions are treated with LES. According to Peng (2005) the hybrid techniques can be subdivided into flow matching and turbulence matching methods. Within the flow matching methods the interface between URANS and LES is explicitly defined. LES filtered equations are solved in the LES region, whereas URANS equations are solved in the URANS domain. The flow parameters (velocities and kinetic energy) are matched at the interface between the URANS and the LES regions. Among the most important contributions to the development of flow matching methods we mention the works of Davidson and Dalstroem (2005), Terracol (2005), Jakirlic et al. (2007), Temmerman et al. (2005) and others. A serious weakness of this approach is the development of robust procedures to set the URANS–LES interface for complicated flow geometries. Within the framework of the turbulence matching method an universal transport equation is solved in the whole computational domain. The stress terms in this equation are treated in different ways in LES and URANS domains. There are various procedures to distinguish between LES and URANS cells. The most popular hybrid method is detached eddy simulation (DES) proposed by Spalart et al. (1997). The original version of this method is based on the classic Smagorinsky LES model and the Spalart–Allmaras (SA) URANS approach. SA is used close to the wall, whereas LES in the rest part of the flow. The switching between the two techniques is smooth and occurs in a ‘‘gray’’ subdomain. There are two major improvements of DES, developed recently. The first one, DDES (delayed DES), has been proposed to detect the boundary layers and to prolong the RANS mode, even if the wall-parallel grid spacing would normally activate the DES limiter (Spalart, 2009). The second one, IDDES (improved DDES), allows one to solve the problems with modeled-stress depletion and log-layer mismatch. For the details see the review Spalart (2009). In spite of a wide application area DES has serious principle limitations thoroughly analyzed by Menter and Egorov (2005). Other versions of the turbulence matching methods using different blending functions to switch the solution between LES and URANS modes were proposed by Peng (2005), Davidson and Billson (2006), Abe and Miyata (2005) and others. A very critical point of the turbulence matching methods is the transition from the time (or ensemble) averaged smooth URANS flow to the oscillating LES flow (see Menter and Egorov, 2005). The oscillations have to appear within a short flow domain in a ‘‘gray zone’’ between LES and URANS. Experience shows that it is extremely difficult to provide a smooth transition of the turbulent kinetic energy passing from the URANS to LES domain. To overcome this problem Schlueter et al. (2002) and Benerafa et al. (2005) used an additional forcing term in the Navier Stokes equation artificially enhancing fluctuations in the gray zone. However, the problem of smooth solution transition from URANS to LES still remains as the main challenge for the turbulence matching methods. This paper presents first results of the development of hybrid methods for ship hydromechanics application undertaken at the

Chair of Modeling and Simulation of the Rostock University. The development proceeds in two directions using both turbulence matching and flow matching techniques. In this paper only a turbulence matching hybrid method is presented with application to prediction of the unsteady wake behind the tanker KVLCC2.

2. Description of the hybrid model Our hybrid model is based on the observation that the basic transport equations have the same form in LES and RANS l

t

@u i @ðu i u j Þ @p n @ðtij þ tij Þ þ ¼ þ @xj @t @xi @xj

ð1Þ

but the interpretation of the overline differs. In LES it means filtering, but in RANS it stands for the Reynolds, or ensemble, averaging. Here we used the standard notation of pn for the pseudo-pressure, and tlij and ttij for the laminar and turbulent stresses respectively. Note that the turbulent stresses are calculated in different ways in LES and URANS regions. The computational domain in our model is dynamically (i.e. at each time step) divided into the LES and URANS regions. A cell of the mesh belongs to one or the other region depending on the relation between the integral length scale L and the extended LES filter D according to the following rule: if L4 D then the cell is in the LES region, if Lo D then the cell is in the URANS region:

ð2Þ

The integral length scale is calculated from the known formula of Kolmogorov and Prandtl with the correction factor 0.168 taken from Schlichting (2000) L ¼ 0:168

k3=2

e

,

ð3Þ

where k is the turbulent kinetic energy and e is the dissipation rate. L varies from one time step to another, which results in varying decomposition of the computational domain into the LES and URANS regions. The extended LES filter is computed as qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð4Þ D ¼ d2max þ d2 , where dmax is the maximal length of the cell edges dmax ¼ maxðdx ,dy ,dz Þ and d ¼ ðthe cell volumeÞ1=3 is the common filter width used in LES. This choice ensures that very flat cells in the boundary layer (for which d  0 but dmax 40) are treated correctly. D depends only on the mesh and it is precomputed only once before the main computation. The LES and URANS regions are shown in Fig. 1. The URANS region is located close to the ship surface and plays the role of a dynamic wall function. In areas of bilge vortices formation, the boundary layer is shedding from the hull and penetrates into the outer flow part. Since the boundary layer is a fine scale flow the procedure (2) recognizes the bilge vortex formation zones as URANS ones. There is a technical issue concerning the cells which

Fig. 1. The division of the computational domain into the URANS (dark) and LES (light) regions at one time instant.

N. Kornev et al. / Ocean Engineering 38 (2011) 1831–1838

are far from the ship hull and where both k and e are small, so large numerical errors are introduced into the integral length scale computed according to Eq. (3). To avoid an irregular distribution of URANS and LES zones, the general rule (2) of the domain decomposition is corrected in such a way that the LES region is switched to URANS one if k is getting less than some threshold. This procedure has no influence on the ship flow parameters since it is used far from the area of the primary interest. We have performed several calculations with different combinations of LES and URANS models to find the most efficient one for the problem under consideration. Among the models we used in our computations are the linear and nonlinear k2e, k2o SST and kev2 f URANS models combined with the simple and dynamic Smagorinsky as well as with the dynamic mixed LES closure models. The experience shows that the most satisfactory results are obtained using the URANS approach based on the kev2 f turbulent model of Durbin (1991) and LES approach based on the Smagorinsky dynamic model. The turbulent stresses ttij are calculated from the Boussinesq approximation using the concept of the turbulent viscosity. The only difference between LES and URANS is the definition of the kinematic viscosity. Within LES it is considered as the subgrid viscosity and calculated according to the dynamic model of Smagorinsky:   1 @uj @ui nSGS ¼ cD d2 9Sij 9, Sij ¼ þ , ð5Þ 2 @xi @xj where Sij is the strain velocity tensor and cD is the dynamic constant. In the URANS region the viscosity is calculated from the turbulent model of Durbin (1991):   k2 nt ¼ min 0:09 ,0:22v2 Tt , ð6Þ

e

2

where v is the wall normal component pffiffiffiffiffiffiffiffi of the stresses and Tt is the turbulent time Tt ¼ maxðk=e,6 n=eÞ.

3. Numerical results for the tanker KVLCC2 The hybrid model described in the previous section has been implemented in the open-source code OpenFOAM based on the finite volume method. The doubled model of the KRISO tanker 1 has KVLCC2 (Kim et al., 2001; Lee et al., 2003 with the scale 58 been chosen for our investigations since it is a well-known benchmark which is widely used in the shipbuilding community (Gothenburg, 2010). The model has length of 5.517 m, breadth of 1 m, draught of 0.359 m and block coefficient of 0.8098. The computations have been carried out on two unstructured 3Dgrids. The coarse grid contains 1.2  106 cells, the fine one—1.8  106 cells. Both grids are coarse in the foreship area (y1þ  10). The grid with 1.2  106 cells has y1þ  226 in the wall region of aftership. The grid with 1.8  106 cells was additionally refined in the aftership area and has y1þ  0:124 in the wall region. This grids were generated by the Ship Model Basin Potsdam (SVA Potsdam) and proved to be an appropriate grid for RANS calculations. The computations on the both grids have been carried out with fixed maximal Courant number of 0.6. The time step was about 0.0008 s for the coarse grid and 0.0005 s for the fine one. For the space discretization, the central difference scheme is used for all terms in the momentum equation whereas the Crank– Nicholson scheme is used for the time discretizations. Steady RANS solutions are used to initialize the flow in the computational domain. The time-averaged solutions has been obtained when the resolved flow reached a statistically steady state (usually it requires the ship way of 3–4 lengths). A typical time

1833

period for statistical averaging takes about 40–50 s, which corresponds to 8–10 lengths of the ship way. Study of the wake has been performed for the constant velocity of 1.047 m/s corresponding to the Reynolds number of Re ¼4.6  106. The Froude number Fn ¼0.142 is small which makes it possible to neglect the water surface deformation effects. 3.1. Estimations of the resolution necessary for a pure ship LES The precise determination of the necessary LES resolution is quite difficult. Estimations presented below are based on the idea that about 80% of the turbulent kinetic energy should be directly resolved and the rest is modeled in a properly resolved LES simulation. Implementation of this idea implies the knowledge of the Kolmogorov Z and the integral length L scales which are used to draw the typical spectra of the full developed turbulence E(k). The wave number kn separating the resolved and modeled turbulence is found from the condition R1 n EðkÞ dk Rk1  0:2: ð7Þ 0 EðkÞ dk The maximum possible cell size is then Dmax ¼ 2p=kn . The scales L and Z are found from the known expression Z ¼ ðn3 =eÞ1=4 and Eq. (3), where the kinetic energy k and the dissipation rate e are taken from RANS simulations using k2e linear model. The ratio l ¼ Dmax =Z is then used as the scale parameter for grid generation. Both lengths vary in space which makes the grid generation procedure very complicated. To roughly estimate the size of the grid we assume that l is constant. We performed different calculations determining l at the two following points: (i) the point where L=Z is maximal in the boundary layer and (ii) the point in the propeller disk where the vorticity x is maximal (region of the concentrated vortex structure). The latter is dictated by the wish to resolve the most intensive vortex flow structures which have the strongest influence on the propeller operation. Since LES application is required in the ship stern area only this part of the computational volume has been meshed. It covers the boundary layer of the stern region starting from the end of the parallel midship section. The thickness of the meshed region has been constant and equal to the maximum boundary layer thickness at the stern dBL . The grid for a pure LES is generated using the following algorithm. The minimum Kolmogorov length Zmin is determined in the near wall region. The cell sizes in x and z directions along the wall are calculated by multiplication of Zmin with the scale parameter l. These sizes remain constant for all cells row in y direction which is normal to the ship surface (see Fig. 2). The cells have at least two equal sizes which is desirable from the point of view of LES accuracy. The choice of the size in y direction is dictated by proper resolution of the boundary layer. Close to the wall this size is chosen from the condition Dw ¼ minðyw , Zmin Þ. Since yw is chosen as the ordinate

Fig. 2. The cell parameters.

1834

N. Kornev et al. / Ocean Engineering 38 (2011) 1831–1838

where y þ ¼1 the first nodes lay deeply in the viscous sublayer. The size in y direction at the upper border of the boundary layer is equal to D1 ¼ lZd , where Zd is the Kolmogorov scale at y ¼ dBL . A simple grading is used in y direction between Dw and D1 . Results of the estimations are as follows: the required grid size ranges from  5 M to  25 M for Re¼ 2.8  106, and from  7 M to  60 M for Re ¼5.8  106. The results vary depending on the value of l in use, so they should be considered as very rough estimations. Together with similar estimations for the nonlinear k2e model these results show that the LES grid should have the order of tens of millions of nodes. Nowadays, the computations with hundred millions and even with a few billions of nodes are becoming available in the research community. However, a numerical study of engineering problems implies usually many computations which have to be performed within a reasonable time with moderate computational resources. In this sense, the results of the present subsection clearly demonstrate that the pure LES is impossible for ship applications so far. To verify that the resolution estimation procedure we used gives meaningful results, it has been applied for turbulent boundary layer (TBL) benchmark. We found from methodical calculations that the pure LES with 1 M cells is quite accurate for prediction of the velocity distribution, TBL thickness, TBL displacement thickness and the wall shear stress. The estimation procedure presented above predicted the necessary resolution around 0.5 M. Therefore, the estimations presented for a ship model are rather lower bound for the resolution required for a pure LES. In the CFD community one can observe tendency to use pure LES without paying any attention to resolution problems. Very often LES is running on typical RANS grids. In fact, such computations can give correct results if the flow structures to be captured are large enough and exist for a long time. In some cases modeling of such structures does not require detailed resolution of boundary layers and a thorough treatment of separation regions. As an example one can mention flows around bluff bodies with predefined separation lines like ship superstructures. Application of underresolved LES for well streamlined hulls should be considered with a great care. First of all, one should not forget that the basic LES subgrid models are derived under the assumption that at least the inertial turbulent subrange is resolved. Second, underresolution of wall region leads to a very inaccurate modeling of the boundary layer, prediction of the separation and overall ship resistance. It is clearly illustrated in the Table 1. The resistance obtained from underresolved LES using the wall function of Werner and Wengle (1991) is less than half of the measured one and that obtained from RANS. Obviously, the application of modern turbulence LES models, more advanced than RANS models, does not improve but even makes the results much worse with the same space resolution. The change from RANS to LES should definitely be followed by the increase of the resolution which results in a drastic increase of the computational costs. These facts underline the necessity of further development towards hybrid methodology. Although in Alin et al. (2010) it has been shown that the accuracy of the resistance prediction using pure LES at a very moderate resolution with y þ  30 can be

improved using special wall functions, the most universal way for the present, to our opinion, is application of hybrid methods.

3.2. Validation Before we start to analyze unsteady effects, we show that our hybrid method predicts averaged flows with the accuracy not worse than that of RANS. Table 1 confirms that the hybrid method works well for ship resistance prediction. Both the overall resistance and the resistance components ratio are in a good agreement with the KRISO measurements (Kim et al., 2001) and RANS. Disagreement between LES results and measurements is due to coarse resolution of the boundary layer. The axial mean velocity field in the propeller plane for the KVLCC2 shown in Fig. 3(a) is compared with the experimental data of KRISO (Lee et al., 2003). The axial velocity ux is normalized to the ship model velocity u0 and the coordinates are normalized to the length between perpendiculars of the ship model Lpp. The mean velocity field is very similar to the experimental one. The lines of the constant velocity have the typical form and reflect the formation of a large longitudinal bilge vortex in the propeller disk. The second longitudinal vortex is formed near the water plane, but it has a much smaller strength compared to the bilge one. More detailed comparison is given in Fig. 4. The axial mean velocity in circumferential direction at the propeller radii r/R¼0.7 and r/R¼ 1 is compared with KRISO experimental data and

Table 1 Results of the resistance prediction using different methods. CR is the resistance coefficient, CP is the pressure resistance and CF is the friction resistance.

KRISO Exp. RANS kev2 f k2o SST-SAS Underresolved LES Hybrid RANS LES

CR

CP (%)

CF (%)

4.11  10  3 4.00  10  3 3.80  10  3 1.70  10  3 4.07  10  3

15 16 18 81 17

85 84 82 19 83

Fig. 3. The mean axial velocity field ux/u0 in the propeller plane computed with different models (right) vs. measurements (left). (a) Hybrid URANS–LES, (b) SA-DDES, and (c) k2o SST  SAS.

N. Kornev et al. / Ocean Engineering 38 (2011) 1831–1838

simulations using the k2o SST-SAS and hybrid models. The agreement in the ranges 0 o y o 401 and 901 o y o 1801 can be considered as quite satisfactory. The results for the range 401 o y o 901 require additional numerical investigations and measurements. However, bearing in mind a big scattering of experimental data in this range, the observed discrepancy is not necessary the sign of modeling weakness. Fig. 5 shows the resolved turbulent kinetic energy normalized to u20 for both computational grids in comparison with experimental data of KRISO. Topologically, the isolines are similar to those of the axial mean velocity shown in the Fig. 3(a). Refinement of the grid leads to a significant improvement of the numerical results. The position of the area with the strongest fluctuations and the magnitudes of these fluctuations are reproduced better when the resolution increases. According to Lee et al. (2003), the uncertainty of the measured TKE is  12%, so our results are quite satisfactory. The OpenFOAM implementation of the detached eddy simulation approach DDES failed to predict the averaged velocity field properly, see Fig. 3(b). Advanced URANS technique k2o SST-SAS provides quite satisfactory results, see Figs. 3(c) and 4.

1835

1 Hz. It corresponds to the Strouhal number of Sh ¼ oB=u0  1 based on the ship model velocity u0 and the ship model breadth B. Similar spectra, calculated with the symmetry condition enforced in the ship symmetry plane (semihull case), revealed no pronounced dominating frequency. Therefore, we ascribe the observed dominating frequency to the interaction between the vortex structures arising at the opposite boards. Like behind bluff bodies this interaction results in the periodical transversal wake oscillations. 3.4. Unsteady effects in the wake

behind the deadwood, attains the maximum at r=R  0:7 and then tends again to zero at large distance from the hull. The frequency analysis allows one to distinguish a dominating frequency around

The most important finding of this work is the observation of the strong unsteadiness of the velocity field in the wake which has not been taken into account in modern engineering methods yet. Fig. 9 shows the history of the axial velocity at the points 3 and 4. The time averaged (mean) velocities are 0.08 and 0.50 at points 3 and 4 respectively. Actually, the velocity changes from  0.06 to 0.23 m/s at the point 3 and from 0.02 to 0.72 m/s at the point 4. Obviously, the estimations of the thrust and cavitation inception are quite different for velocities 0.50 and, say, 0.02 m/s. In the case of ux ¼0.02 the thrust is much larger and the cavitation is more probable. The variation range of the axial velocity is characterized by the width of the probability density function (PDF) of ux (see Fig. 10). It should be noted that the PDF distributions are generally asymmetric. The skewness of the PDF is sufficiently negative at points 2, 4 and 5. It means that the PDF at ux þa decays faster than that at ux  a when a increases. This results in long tails at small velocities which means that very small velocities are more probable than very large ones. Therefore, the temporal variation of the thrust should be very asymmetric with respect to its mean value. Sufficiently large thrust values and cavitation inception corresponding to small velocities are very probable. We emphasize that the strong unsteadiness of the wake has been obtained only using the present hybrid method. URANS and SAS show good performance for the time averaged values. Unsteadiness predicted by these methods is either unrealistically weak (SAS) or

Fig. 4. Circumferential distribution of the mean axial velocity field in the propeller plane. W—k2o SST  SAS, J—DSMþ V2F, solid line—KRISO experiments for the specified r/R. (a) r/R ¼0.7 and (b) r/R¼ 1.0.

Fig. 6. Positions of probe points. R is the propeller radius.

3.3. Spectra of the axial velocity The one-dimensional spectrum of the turbulent kinetic energy (TKE) E11 is presented in Fig. 7 at the set of points shown in Fig. 6. All of the presented spectra inside of the propeller disk show the slope of 5=3 confirming the fact that the inertial subrange of the TKE has been properly resolved even at relatively coarse resolution typical rather for RANS than for LES. Obviously, the vortex structures in the propeller disk are so large that they can be correctly reproduced on the given grid. The kinetic energy qffiffiffiffiffiffiffi represented by the standard deviation s ¼ u02 x is close to zero

Fig. 5. Resolved turbulent kinetic energy k ¼ r=2ðu0x u0x þ u0y u0y þ u0z u0z Þ=u20 multiplied with 103 in the propeller plane. Numerics (right-half of each figure) versus measurement (left-half). (a) Computational grid 1.2  106 cells and (b) computational grid 1.8  106 cells.

1836

N. Kornev et al. / Ocean Engineering 38 (2011) 1831–1838

1

2

3

4

5

6

Ux (m/s)

Fig. 7. One-dimensional spectra E11 at different points along the line y ¼ 901. Solid line E11  o5=3 . Number of samples ¼ 69,150, the time of statistical analysis is treal ¼ 53 s with Dt ¼ 7:7  104 s. Cut-off frequency is  4 kHz.

1

6 5 4 2 3 1

0.6 0.2

-0.2 0

5

10

15

t (s) Fig. 8. History of ux obtained using k2o SST-SAS model.

even absent at all (various URANS models). It can be seen in Fig. 8 where the velocity history obtained using k2o SST-SAS is presented. There are no oscillations within 15 s of the simulation.

4. Estimation of the thrust variations due to unsteady wake flow To roughly estimate the impact of flow unsteadinesses described above we use a simple profile theory neglecting the induced velocities. The thrust coefficient KT on a blade is given by the integral Z 1 1 dK T KT ¼ dr, ð8Þ 4 r hub dr where the function under the integral reads as follows (see, for instance, Artjushkov et al., 1984): dK T b  vR 2 ¼ CL cos bð1e tan bÞ, ð9Þ D n2R dr with r being the radial coordinate, r ¼ r=R, CL the lift coefficient of the profile, b the local profile chord, vR the incident velocity on profile calculated neglecting the induced components, i.e. qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi vR ¼ ðOrÞ2 þ v20 , O is the angular velocity of the propeller, e is the reciprocal lift-to-drag ratio and b is the angle between vR and the velocity component Or. The coefficients CL and e are calculated using the approximations proposed by Mishkevich (see Artjushkov et al., 1984) CL ¼ 2pmða þ2wf Þ,

e¼

0:05808ð1 þ 2:3tÞ CL Re1=7

,

ð10Þ

where m and n are given by the following expressions:

m ¼ ð1 þ 0:87tÞð1e0:0691 þ 12:46t 0:1855 ln Re Þ, w ¼ 1:015 1 þ

tðt0:05Þ ð0:04664 ln Re0:4378Þ2

! ð11Þ

and the Reynolds number is defined in the usual way as Re ¼ bvR =n. Here t denotes the profile thickness, t ¼ t=b, f the curvature f ¼ f =b, a the angle of attack a ¼ jb and j the propeller pitch angle. The propeller radius R¼ 0.085 m has been taken from Kim et al. (2001). The chosen rotation frequency n ¼5.9 Hz corresponds to the optimal advance ratio. The radial distributions of the thickness t, the curvature f, the profile chord b and the pitch angle j are chosen the same as for the propeller DTRC 4119 without skew. The thrust coefficient of the whole propeller is the sum of the contributions of all its blades. To perform the thrust calculations the velocity field has been precomputed within 53 s, stored and then used for the calculation of the thrust. In Figs. 11–13 we presented the results for dK T =dr calculated at r/R¼0.7 with consideration of the wake unsteadinesses and compared them with the case of the frozen wake calculated by averaging the unsteady wake. It is seen that the standard deviation of dK T =dr decays when the blade number increases, which is in agreement with the present experience. In the case of the frozen wake the decay is faster than in the case of the unsteady wake. It is noticeable that the thrust deviation for the number of blades Z 4 4 increases substantially when unsteady effects are considered (compare the solid and dashed lines). The spectral analysis of the slope dK T =dr performed for a five-blade propeller revealed the presence of the dominating frequencies which are proportional to nZ. We see again that this result is in a well agreement with the previous observations. The second mode, which has the frequency 2nZ (59 Hz in our case), is not readable in the frozen wake calculations since high frequencies are filtered out by averaging the original unsteady wake. An interesting peculiarity of the unsteady wake result (Fig. 12(a)) is the presence of an additional low-frequency peak at o  0:7 Hz which is, obviously, the consequence of the dominating oscillations in the wake with the frequency of about 1 Hz. The most important result is presented in Figs. 13 and 14 which shows that the unsteady flow causes stronger variations of the thrust than those for the frozen wake. In the frozen wake the oscillations are regular and strictly periodic, whereas they are only quasi-periodic and irregular in the unsteady wake case. The averaged thrust is in the unsteady case approximately 2.5% larger than in the case of the frozen wake although the averaged wake velocities are exactly the same. This is the effect of the oscillations and the nonlinearity caused by the velocity squared v2R in Eq. (9). The variation of the full thrust calculated from Eq. (8) and related to the mean thrust DKT =KT0 is 10% for the frozen wake and 26% for the unsteady wake. The PDF distribution of the thrust (see Fig. 14)

N. Kornev et al. / Ocean Engineering 38 (2011) 1831–1838

1837

Fig. 9. History of ux at points 3 and 4. Dotted lines are time averaged values.

Ux = 0

1

10

30

8

20

6

Ux = 0.27

2

0

0

-0.08

-0.04

0

0.04

Ux = 0.51

6

0.1 4

4 2 0

20

0.2

0.3

0.4 5

Ux = 0.84

0.5

0 80

15

60

10

40

5

20

0.7

0.05

0.1

0.15

0.2 6

Ux = 0.953

0

0 0.3

3

3

2

0

Ux = 0.09

6

4

10

9

0.78

0.82

0.86

0.9

0.94

0.95

0.96

0.97

Fig. 10. The probability density function of Ux at different points. Solid lines are corresponding Gaussian pdfs. Number of samples ¼ 69,150, treal ¼ 53 s, Dt ¼ 7:7  104 s.

Standard deviation

0.05 0.04 0.03 0.02 0.01 3

4

5 6 Number of blades

7

8

Fig. 11. The standard deviation of dK T =dr ðr ¼ 0:7Þ as a function of the blade number for the unsteady (solid line) and averaged (frozen) wakes (dashed line).

m =1

1.5

m =2

m =1

1.5

1

1

0.5

0.5

m =2

0

0 0

10 20 30 40 50 60 ω

0

10 20 30 40 50 60 ω

Fig. 12. The spectra of dK T =dr ðr ¼ 0:7Þ. Vertical lines are the frequencies proportional to nZ. (a) Unsteady wake and (b) averaged (frozen) wake.

shows that the thrust values slightly larger than the mean value are more frequent than slightly smaller ones.

5. Conclusion Unsteady loadings are of a great interest for many practical applications in marine engineering. Modeling of these

phenomena requires new theoretical tools capable of resolving concentrated vortex structures which are responsible for large scale flow fluctuations. The focus of the present paper is the study of unsteady effects in the wake of ships with large block coefficients. The estimations presented in this paper show that the pure LES demands a huge resolution which makes the LES application for shipbuilding purposes impractical even at small Re numbers typical for ship models. So far, the only reasonable solution is the application of hybrid methods. In this paper we presented a hybrid method based on the combination of the dynamic Smagorinsky SGS model (DSM) and the kev2 f URANS approach. The method is applied to the calculation of the resistance and the wake flow of the tanker KVLCC2. The hybrid method provides very good results for the resistance. Also the fields of the mean axial velocity and the turbulent kinetic energy agree well with the measurement in the propeller plane. The hybrid method predicts the unsteadiness of the wake flow. On the contrary, all URANS models including SAS as well as DES approach proved to be not capable of reproducing the unsteady character of the wake flow in the propeller disk at least on the relatively coarse grids used in this paper. Analysis shows that the instantaneous velocities deviate sufficiently from the mean values which are usually used as the estimated velocities in modern engineering methodologies. This fact can negatively influence the accuracy of propulsion and unsteady loads prediction. The probability density function of the axial velocity is sufficiently asymmetric at points within the propeller disk. As a result, the temporal variation of the thrust is very asymmetric with respect to its mean value. Thrust variation in unsteady wake is more than twice as large as that in the time averaged (frozen) wake. The results of the present study point out that the unsteadiness in the wake behind full ships can be very large and should be taken into account when propulsion and unsteady loadings are determined.

1838

N. Kornev et al. / Ocean Engineering 38 (2011) 1831–1838

Fig. 13. Time dependency of dK T =dr ðr ¼ 0:7Þ for the unsteady wake (solid line) and for the frozen wake (dashed line). Z¼ 5.

25 20 15 10 5 0 0.38

0.40

0.42

0.44

0.46

0.48

KT Fig. 14. Probability density function of the whole thrust. The vertical line shows the mean value. Z¼ 5.

Acknowledgment The authors acknowledge gratefully the support of the German Federal Ministry of Economics and Technology (BMWi). This work has been performed within the framework of the project ShipLES monitored by Dr. A. Nitz. LES calculations have been performed on IBM pSeries 690 Supercomputer at the North German Alliance for the Advancement of High-Performance Computing (HLRN). The authors wish to acknowledge the community developing the OpenFoam.

References Abe, K., Miyata, Y., 2005. An investigation of hybrid LES/RANS models for predicting flow fields with separation. In: Proceedings of the 4th International Symposium on Turbulence and Shear Flow Phenomena, pp. 1153–1158. Alin, N., Bensow, R.E., Fureby, C., Huuva, T., Svennberg, U., 2010. Current capabilities of DES and LES for submarines at straight course. J. Ship Res. 54, 184–196. Artjushkov, L.S., Achkinadse, A.S., Russetsky, A.A., 1984. Marine Propellers. Sudostroenie (in Russian).

Benerafa, Y., Ducros, F., Sagaut, P., 2005. RANS/LES coupling using a forcing term approach. In: Proceedings of the 4th International Symposium on Turbulence and Shear Flow Phenomena, pp. 1141–1146. Davidson, L., Billson, M., 2006. Hybrid LES–RANS using synthesized turbulent fluctuations for forcing in the interface region. Int. J. Heat Fluid 27, 1028–1042. Davidson, L., Dalstroem, S., 2005. Hybrid RANS-LES: an approach to make LES applicable at high Reynolds number. Int. J. Comput. Fluid Dynamics 19 (6), 415–427. Durbin, P.A., 1991. Near-wall turbulence closure modeling without ‘‘damping functions’’. Theor. Comput. Fluid Dyn. 3, 1–13. Gothenburg, 2010. A Workshop on CFD in Ship Hydrodynamics. Chalmers University of Technology, Gothenburg. Jakirlic, S., Saric, S., Kniesner, B., Kadavelil, G., Basara, B., Chaouat, B., 2007. SGS modelling in LES of wall-bounded flows using transport RANS models: from a zonal to a seamless hybrid LES/RANS method. In: Proceedings of the 6th International Symposium on Turbulence and Shear Flow Phenomena, pp. 1057–1062. Kim, S., Rhee, S.H., 2002. Assessment of eight turbulence models for a threedimensional boundary layer involving crossflow and streamwise vortices. Fluent Technical Notes 165, 1–25. Kim, W.J., Van, S.H., Kim, D.H., 2001. Measurement of flows around modern commercial ship models. Exp. Fluids 31, 567–578. Kornev, N., 2010. Propeller Theory. Shaker Verlag. Lee, S.J., Kim, H.R., Kim, W.J., Van, S.H., 2003. Wind tunnel tests on flow characteristics of the KRISO 3600 TEU containership and 300 KVLCC doubledeck ship models. J. Ship Res. 47, 24–38. Menter, F., Egorov, J., 2005. Turbulence models based on the length-scale equation. In: Proceedings of the 4th International Symposium on Turbulence and Shear Flow Phenomena, pp. 941–946. NuTTS, 2010. 13th Numerical Towing Tank Symposium, Duisburg. Peng, S., 2005. Hybrid RANS–LES modeling based on zero- and one-equation models for turbulent flow simulation. In: Proceedings of the 4th International Symposium on Turbulence and Shear Flow Phenomena, pp. 1159–1164. Schlichting, H., 2000. Boundary Layer Theory. Springer. Schlueter, J.U., Pitsch, H., Moin, P., 2002. Consistent boundary conditions for integrated LES/RANS simulations: LES outflow conditions. AIAA J. 3121, 1–8. Spalart, P.R., Jou, W.H., Strelets, M., Allmaras, S.R., 1997. Comments on the feasibility of LES for wings and on a hybrid RANS/LES approach. In: Advances in LES/DNS: Proceedings of the First AFOSR International Conference on DNS/ LES, pp. 137–148. Spalart, P.R., 2009. Detached-Eddy simulation. Annu. Rev. Fluid Mech. 41, 181–202. Temmerman, L., Hadziabdec, M., Leschziner, M., Hanjalic, K., 2005. A hybrid twolayer URANS–LES approach for large-eddy simulation at high Reynolds numbers. Int. J. Heat Fluid 26, 173–190. Terracol, M., 2005. Airframe noise prediction by mean of a zonal RANS/LES approach. In: Proceedings of the 4th International Symposium on Turbulence and Shear Flow Phenomena, pp. 1165–1169. Werner, H., Wengle, H., 1991. Large-eddy simulation of turbulent ow over and around a cube in a plate channel. In: 8th Symposium on Turbulent Shear Flows.