Advanced
Investigation on the wall function implementation for the prediction of ship resistance
Investigation on the wall function implementation for the prediction of ship resistance
International Journal of Naval Architecture and Ocean Engineering. 2013. Mar, 5(1): 33-46
Copyright ©2013, The Society of Naval Architects of Korea
  • Published : March 31, 2013
Download
PDF
e-PUB
PubReader
PPT
Export by style
Share
Article
Author
Metrics
Cited by
TagCloud
About the Authors
Sunho Park
Dept. of Naval Architecture and Ocean Engineering, Seoul National University, Seoul, Korea
Se Wan Park
Dept. of Naval Architecture and Ocean Engineering, Seoul National University, Seoul, Korea
Shin Hyung Rhee
Professor, Dept. of Naval Architecture and Ocean Engineering, Research Institute of Marine Systems Engineering, Seoul National University, Seoul, Korea
shr@snu.ac.kr
Sang Bong Lee
Maritime Research Institute, Hyundai Heavy Industries Co., LTD., Ulsan, Korea
Jung-Eun Choi
Maritime Research Institute, Hyundai Heavy Industries Co., LTD., Ulsan, Korea
Seon Hyung Kang
Maritime Research Institute, Hyundai Heavy Industries Co., LTD., Ulsan, Korea
Abstract
A computational fluid dynamics (CFD) code, dubbed SNUFOAM, was developed to predict the performance of ship resistance using a CFD tool kit with open source libraries. SNUFOAM is based on a pressure-based cellcentered finite volume method and includes a turbulence model with wall functions. The mesh sensitivity, such as the skewness and aspect ratio, was evaluated for the convergence. Two wall functions were tested to solve the turbulent flow around a ship, and the one without the assumption of the equilibrium state between turbulent production and dissipation in the log law layer was selected. The turbulent flow around a ship simulated using SNUFOAM was compared to that by a commercial CFD code, FLUENT. SNUFOAM showed the nearly same results as FLUENT and proved to be an alternative to commercial CFD codes for the prediction of ship resistance performance.
Keywords
NOMENCLATURE
  • GkProduction of turbulent kinetic energy [kg/ms3]
  • IUnit tensor [-]
  • kTurbulent kinetic energy [m2/s2]
  • PPressure [Pa]
  • pEstimated order of accuracy of the computational method [-]
  • RERichardson extrapolated value [-]
  • ReReynolds number [-]
  • uτFriction velocity [m/s]
  • U∞Free stream velocity [m/s]
PPT Slide
Lager Image
  • Velocity vector [m/s]
  • ypDistance from wall to cell center [m]
  • εTurbulent dissipation rate [m2/s3]
  • μViscosity [kg/m-s]
  • μeffEffective viscosity [kg/m-s]
  • μtTurbulent viscosity [kg/m-s]
  • νtTurbulent eddy viscosity [m2/s]
  • ρDensity [kg/m3]
  • κSurface roughness [-]
PPT Slide
Lager Image
  • Turbulent stress tensor [Pa]
  • τwWall shear stress [Pa]
  • φComputational solution variable for uncertainty assessment [-]
INTRODUCTION
With the rapid development of CFD techniques, computational analyses complement experiments in many parts of the hullform design. Korea Ocean Research and Development Institute (KORDI) developed a CFD code for the prediction of ship resistance and propulsion performances (Kim et al., 2002). Many shipbuilding companies in Korea adopted the code as their CFD tool for hull form development. FLOWTECH also developed a CFD code for ship resistance performance prediction (Leer-Andersen and Larsson, 2003). General purpose CFD codes, such as FLUENT, STAR-CCM+, and CFX, are also used in shipbuilding companies, research institutes, and universities (Rhee et al., 2005; Rhee, 2009; Park and Chun, 2009; Choi et al., 2010; Seo et al., 2010, Park et al., 2010). Specialized CFD codes for hull form development, such as the ones developed by KORDI and FLOWTECH, are equipped with limited computational methods and their applications are limited to common ship hull shapes. On the other hand, general purpose CFD codes can handle non-conventional ships and provide various computational methods. However, their purchase and maintenance prices are quite high. Especially, annual licensing fee for parallel computing is another burden. Due to these reasons, interest in free CFD codes has increased recently. Among free CFD codes, Open FOAM, which is an object-oriented, open source CFD tool kit draws huge interest recently (Jasak, 2009).
Many CFD codes and libraries have been developed using the OpenFOAM platform. Favero et al. (2009) developed a new tool for CFD simulation of viscoelastic fluids. Kunkelmann and Stephan (2009) implemented and validated a nucleate boiling model using OpenFOAM platform. Bensow and Bark (2010) employed an approach to simulate dynamic cavitation behavior based on the large eddy simulation, using an implicit approach for the subgrid terms and applied to a propeller flow. Kissling et al. (2010) developed a coupled pressure based solution algorithm to model interface capturing problems for multiphase flows. Silva and Lage (2011) coupled a multiphase flow formulation with the population balance equation solution by the direct quadrature method of moments. Dittmar and Ehrhard (2011) developed a phase change model library including the contact line evaporation model and on the conjugate heat transfer between solid and fluid. Petit et al. (2011) investigated the swirling flow with helical vortex breakdown in a conical diffuser. Park and Rhee (2012, 2013) developed a CFD code (SNUFOAM-Cavitation) and applied it to super- and cloud cavitations. Many researchers have developed new libraries and CFD codes for their purpose. However, studies on free CFD codes in ship hydrodynamics are yet to be reported, even though the dependency for commercial CFD codes have been intensified over many years.
Therefore, in the present study, a Reynolds-averaged Navier-Stokes (RANS) equations solver that couples the velocity and pressure was developed based on a cell-centered finite volume method. The developed CFD code, termed SNUFOAM, was based on the OpenFOAM (Jasak, 2009), and the authors implemented wall function (Launder and Spalding, 1972) library.
The present study focused on the turbulent flow around a surface ship. The objectives were therefore (1) to develop and validate SNUFOAM, and (2) to understand the turbulent flow around a surface ship by comparing the results with a comercial CFD code.
The paper is organized as follows. The description of the physical problem is presented first, and this is followed by the computational methods. The computational results are then presented and discussed. Finally, a summary and conclusions are provided.
COMPUTATIONAL METHODS
- Mathematical modeling
The equations for the mass and momentum conservation were solved to obtain the velocity and pressure fields. The equation for the conservation of mass, or the continuity equation, can be written as
PPT Slide
Lager Image
where
PPT Slide
Lager Image
is the velocity vector.
The equation for the conservation of momentum can be written as
PPT Slide
Lager Image
where P is the static pressure and the turbulent stress tensor (
PPT Slide
Lager Image
), is given by
PPT Slide
Lager Image
with the second term on the right-hand side representing the volume dilation effect, where μeff = μ + μt , μ is the viscosity, and the subscripts eff and t denote “effective” and “turbulent,” respectively. I is the unit tensor.
Once the Reynolds averaging approach for turbulence modeling is applied, the unknown term, i.e., the Reynolds stress term, is related to the mean velocity gradients by the Boussinesq hypothesis, as follows
PPT Slide
Lager Image
The standard k-ε turbulence model, which is based on the Boussinesq hypothesis with transport equations for the turbulent kinetic energy ( k ) and its dissipation rate ( ε ) was adopted for turbulence closure (Shih et al., 1995). The turbulent viscosity ( μt ) was computed by combining k and ε via ρCμk2 , and the turbulent kinetic energy and its rate of dissipation are obtained from the transport equations, which are
PPT Slide
Lager Image
PPT Slide
Lager Image
where Cμ was an empirical constant of 0.09. Here, the model constants C , C2 , and C , were 1.44, 1.92, and 1.0, respectively. C1 was
PPT Slide
Lager Image
The turbulent viscosity was used to calculate the Reynolds stresses to close the momentum equations. The wall function was used for the near-wall treatment (Pope, 2000; Launder and Spalding, 1972). A detailed description of the wall functions is as followed.
From the assumption of the equilibrium state ( Gk = ε ), the local solution for ε at the wall yields,
PPT Slide
Lager Image
where u τ is the friction velocity, κ is the surface roughness, and yp is the distance from the wall surface to cell center. The eddy viscosity at the wall is defined as
PPT Slide
Lager Image
and
PPT Slide
Lager Image
which leads to the expression for u τ in terms of k .
PPT Slide
Lager Image
From the definition of
PPT Slide
Lager Image
the formulation of the production of turbulent kinetic energy in the turbulent kinetic energy equation, proposed by Pope (2000), is expressed, as follows
PPT Slide
Lager Image
where τw is the wall shear stress. The CFD code including the equation (11) was named SNUFOAM-WF1.
The wall shear stress of the log-law is commonly obtained (Launder and Spalding, 1972; Craft et al., 2002).
PPT Slide
Lager Image
Thus,
PPT Slide
Lager Image
From the definition of
PPT Slide
Lager Image
the formulation of the production of turbulent kinetic energy in the turbulent kinetic energy equation, proposed by Launder and Spalding (1972), is expressed, as follows
PPT Slide
Lager Image
The CFD code including the equation (14) was named SNUFOAM-WF2.
- Numerical methods
A pressure-based cell-centered finite volume method was employed along with a linear reconstruction scheme that allows the use of computational cells of arbitrary shapes. The solution gradients at the cell centers were evaluated by the least-square method. The convection terms were discretized using the van Leer scheme (van Leer, 1979), and for the diffusion terms, a central differencing scheme was used. The velocity-pressure coupling and overall solution procedure were based on a pressure-implicit with splitting order (PISO) type segregated algorithm (Issa, 1985) adapted to an unstructured grid. The discretized algebraic equations were solved using a pointwise Gauss-Seidel iterative algorithm, while an algebraic multi-grid method was employed to accelerate solution convergence.
RESULTS AND DISCUSSION
- Problem description
KRISO container ship (KCS), a 3,600 TEU container carrier, was selected as the object ship. The principal particulars are described in Table 1 and the body plan is shown in Fig. 1 (Kim et al., 2001).
Principal particulars of KCS.
PPT Slide
Lager Image
Principal particulars of KCS.
PPT Slide
Lager Image
Body plan of KCS.
In the Cartesian coordinate system adopted here, the positive x -axis was in the streamwise direction, the positive y -axis was in the starboard direction, and the positive z -axis was in the upward direction. Only a half of the ship was modeled due to the y -axis symmetry. The solution domain extent shown in Fig. 2 was -0.5 ≤ x/l ≤ 1, 0 ≤ y/l ≤ 1, and -1 ≤ z/l ≤ 0. Here, l indicates the length of the ship. The upstream inlet and far-field boundary was specified as the Dirichlet boundary condition, i.e., with a fixed value of the velocity. On the downstream exit boundary, the reference pressure with the extrapolated velocity was applied. The free-surface was not considered, i.e. a double body model without the free-surface was considered, and thus, a slip condition was applied on the top boundary. No-slip condition was applied on the ship surface. A symmetric condition was applied on the center plane boundary.
PPT Slide
Lager Image
Boundary conditions.
The mesh generation was carried out using Gridgen V15.17, commercially available software. A single-block mesh of 370,000 hexahedral cells was used as shown in Fig. 3 . Along the length of the ship, 180 cells were used, while along the girth of the ship, 45 cells were applied.
PPT Slide
Lager Image
Solution meshes.
- Uncertainty assessment
To evaluate the numerical uncertainty in the computational results, the concept of the grid convergence index (GCI) was adopted. Three levels of mesh resolution were considered for the solution convergence of the resistance coefficient (C T ) and pressure coefficients at the fore perpendicular (FP) and mid ship with the design water-plane level.
The order of accuracy can be estimated as
PPT Slide
Lager Image
where ϕcoarse , ϕmedium , and ϕfine are solutions at the coarse, medium, and fine levels, respectively. The Richardson extrapolated value (RE) and the convergence index (CI) were also calculated by Equations (16) and (17), respectively.
PPT Slide
Lager Image
PPT Slide
Lager Image
where
PPT Slide
Lager Image
and where r is the effective mesh refinement ratio of
PPT Slide
Lager Image
with the cell count, N , and the number of dimensions, D . Table 2 summarizes the numerical uncertainty assessment results. Overall, the solutions show good mesh convergence behavior with errors from the corresponding RE of less than 0.6%. Table 2
Numerical uncertainty assessment.
PPT Slide
Lager Image
Numerical uncertainty assessment.
- Validation
Steady computations were done for the Reynolds number (based on U of 2.196 m/s and the length of the ship of 7.2786 m ) of 1.4 × 10 7 and the Froude number (based on U and the length of the ship) of 0.26. The residual history of the computation shown in Fig. 4 indicates that the computation is diverged. Diverged results were investigated to know the source of the divergence. Fig. 5 shows the pressure coefficient ( CP = ( P−Po )/0.5 ρU 2 ) distribution around the inlet and far-field boundaries. The pressure fluctuation was observed in high aspect ratio cells and caused the divergence because the computation was an elliptic type problem. The solution divergence was caused by a floating-point error, which led to a mass imbalance. A numerical method using an averaged value of neighboring cells was utilized to prevent the divergence. In this paper, a new mesh was selected instead of a numerical method to converge the solution. In general, high aspect ratio cells were used in external flows, such as ship hydrodynamics and external aerodynamics. Thus, special care is required when a mesh for external flows is planned. A new mesh was generated considering the cell characteristics, as shown in Fig. 6 . The maximum level of aspect ratio decreased from 750 to 160, and nondimensionalized maximum skewness level decreased 0.96 to 0.93. The residual history with the new mesh is shown in Fig. 7 . The results showed good convergence.
PPT Slide
Lager Image
Residual history.
PPT Slide
Lager Image
Computation errors at high aspect ratio cells.
PPT Slide
Lager Image
Changed solution meshes around the bow (maximum aspect ratio: 160).
PPT Slide
Lager Image
Residual history with changed meshes.
Firstly, computations were done with the new mesh using SNUFOAM-WF1. The validation of SNUFOAM was conducted by comparing with the results of a commercial CFD code, FLUENT, which showed good results in ship hydrodynamics (Rhee et al., 2005; Rhee, 2009; Park and Chun, 2009; Choi et al., 2010; Seo et al., 2010). Fig. 8 shows the pressure coefficient, nondimensionalized turbulent kinetic energy ( k /U 2 ) and nondimensionalized turbulent dissipation rate ( ε L /U 3 ) contours around the bow. The pressure coefficient contours were the similar in the results of both CFD codes, while the nondimensionalized turbulent kinetic energy and nondimensionalized turbulent dissipation rate contours were different. The nondimensionalized turbulent kinetic energy and nondimensionalized turbulent dissipation rate were maintained zero around the bow in SNUFOAMWF1. It was observed that the results computed using SNUFOAM-WF1 needed a greater entrance length to develop the turbulence, which was believed to be influenced by the wall function. Generally, the turbulent kinetic energy was larger than the turbulent dissipation rate in the bow region and the turbulent dissipation was larger than the turbulent production in the stern region. The assumption of the equilibrium state between the turbulent production and dissipation in the log law layer (Pope, 2000) delayed the turbulent production in the bow region.
PPT Slide
Lager Image
Comparison of SNUFOAM-WF1 (left) and FLUENT (right) around the bow.
The same flow was also simulated using SNUFOAM-WF2, which did not include the assumption of the equilibrium state. The equation (14) was calculated by substituting the friction velocity
PPT Slide
Lager Image
which was derived from the assumption of the equilibrium state, to the equation (11). Fig. 9 shows the pressure coefficient, nondimensionalized turbulent kinetic energy, and nondimensionalized turbulent dissipation rate contours around the bow. The nondimensionalized turbulent kinetic energy and nondimensionalized turbulent dissipation rate computed using SNUFOAM-WF2 showed the nearly same development with those computed using FLUENT. The wall function (Launder and Spalding, 1972), which did not include the assumption of the equilibrium state, improved the turbulence development around the bow. Thus the turbulent production was larger than the turbulent dissipation in the bow region. Another reason could be numerical errors. The velocity gradient, which was related with the wall shear stress
PPT Slide
Lager Image
plays an important role here. Equation (14) includes the velocity gradient twice, while Equation (11) includes the velocity gradient only once.
PPT Slide
Lager Image
Comparison of SNUFOAM-WF2 (left) and FLUENT (right) around the bow.
Fig. 10 shows the pressure coefficient, nondimensionalized turbulent kinetic energy, and nondimensionalized turbulent dissipation rate contours around the stern. The pressure coefficient, nondimensionalized turbulent kinetic energy, and nondimensionalized turbulent dissipation rate contours computed using SNUFOAM-WF2 showed quite close to the results computed using FLUENT.
PPT Slide
Lager Image
Comparison of SNUFOAM-WF2 (left) and FLUENT (right) around the stern.
CONCLUDING REMARKS
A CFD code was developed using the OpenFOAM (Jasak, 2009), and the authors implemented wall function (Launder and Spalding, 1972) library. The mesh sensitivity was evaluated for the skewness and aspect ratio of the cells. SNUFOAM showed good convergence with low aspect ratio cells. The results computed using SNUFOAM-WF1 needed the long entrance length for the development of the turbulent boundary layer in the bow, and was different from those computed using FLUENT. SNUFOAM- WF2 was developed including another wall function, which did not include the assumption of the equilibrium state between the turbulent kinetic energy and turbulent dissipation rate in the log law layer. The results computed using SNUFOAMWF2 showed quite close to those computed using FLUENT. It is expected that SNUFOAM including the volume fraction transport equation (Park and Rhee, 2011) can substitute a commercial CFD code for the prediction of ship resistance performance in ship hydrodynamics.
Acknowledgements
This work was supported by the World Class University Project (R32-2008-000-10161-0), Multi-Phenomena CFD Research Center (20090093129) and the National Research Foundation of Korea (2010-0022835) funded by the Ministry of Education, Science and Technology of the Korea government.
References
Bensow R.E. , Bark G. 2010 Implicit LES predictions of the cavitating flow on a propeller. Journal of Fluids Engineering 132 (4) 041302.1 - 10
Choi J.-E. , Min K.-S. , Kim J.H. , Lee S.B. , Seo H.W. 2010 Resistance and propulsion characteristics of various commercial ships based on CFD results. Ocean Engineering 37 (7) 549 - 566
Craft T.J. , Gerasimov A.V. , Iacovides H. , Launder B.E. 2002 Progress in the generalization of wall-function treatments. International Journal of Heat and Fluid Flow 23 (2) 148 - 160
Dittmar I. , Ehrhard P. 2011 Numerical study of liquid/liquid slug flow in a capillary microreactor. Proceedings in Applied Mathematics and Mechanics 11 (1) 617 - 618
Favero J.L. , Secchi A.R. , Cardozo N.S.M. , Jasak H. 2009 Viscoelastic flow simulation: Development of a methodlogy of analysis using the software OpenFOAM and differential constitutive equations. Computer Aided Chemical Engineering 27 915 - 920
Issa R,I. 985 Solution of implicitly discretized fluid flow equations by operator splitting. Journal of Computational Physics 62 40 - 65
Jasak H. 2009 OpenFOAM: Open source CFD in research and industry. International Journal of Naval Architecture and Ocean Engineering 1 (2) 89 - 94
Kim W.J. , Van S.H. , Kim D.H. 2001 Measurement of flows around modern commercial ship models. Experiments in Fluids 31 (5) 567 - 578
Kim W.J. , Kim D.H. , Van S.H. 2002 Computational study on turbulent flows around modern tanker hull forms. International Journal for Numerical Methods in Fluids 38 (4) 377 - 406
Kissling K. , Springer J. , Jasak H. , Schutz S. , Urban K. , Piesche M. 2010 A coupled pressure based solution algorithm based on the volume-of-fluid approach for two or more immiscible fluids. 5th European Conference on Computational Fluid Dynamics Lisbon, Portugal 14-17 June 2010
Kunkelmann C. , Stephan P. 2009 CFD simulation of boiling flows using the volume-of fluid method within Open FOAM. Numerical Heat Transfer, Part A: Applications 56 (8) 631 - 646
Launder B.E. , Spalding D.B. 1972 Lectures in mathematical models of turbulence. Academic Press UK
Leer-Andersen M. , Larsson L. 2003 An experimental/numerical approach for evaluating skin friction on full-scale ships with surface roughness. Journal of Maritime Science and Technology 8 (1) 26 - 36
Park D.-W. , Chun H.-H. 2009 Design practice for the stern hull form of a twin-skeg ship. Journal of Maritime Science and Technology 14 (3) 310 - 321
Park S. , Heo J. , Yu B.S. 2010 Numerical study on the gap flow of a semi-spade rudder to reduce gap cavitation. Journal of Marine Science and Technology 15 (1) 78 - 86
Park S. , Rhee S.H. 2011 CFD code development using open source libraries for shipbuilding and marine engineering industries. Journal of the Society of Naval Architects of Korea 49 (2) 151 - 157
Park S. , Rhee S.H. 2012 Computational analysis of turbulent super-cavitating flow around a two-dimensional wedgeshaped cavitator geometry. Computers & Fluids 70 73 - 85
Park S. , Rhee S.H. 2013 Numerical analysis of the three-dimensional cloud cavitating flow around a twisted hydrofoil. Fluid Dynamics Research 45 (1) 015502.1 - 01550220
Petit O. , Bosioc A.I. , Nilsson H. 2011 Unsteady simulations of the flow in a swirl generator, using OpenFOAM. International Journal of Fluid Machinery and Systems 4 (1) 199 - 208
Pope S.B. 2000 Turbulent flows. Cambridge University Press UK
Rhee S.H. , Kawamura T. , Li H. 2005 Propeller cavitation study using an unstructured grid based Navier-Stoker Solver. Journal of Fluids Engineering 127 (5) 986 - 994
Rhee S.H. 2009 Unsteady Reynolds averaged Navier-stokes method for free-surface wave flows around surface-piercing cylindrical structures. Journal of Waterway, Port, Coastal, and Ocean Engineering 135 (4) 135 - 143
Seo J.H. , Seol D.M. , Lee J.H. , Rhee S.H. 2010 Flexible CFD meshing strategy for prediction of ship resistance and propulsion performance. International Journal of Naval Architecture and Ocean Engineering 2 (3) 139 - 145
Shih T.H. , Liou W.W. , Shabbir A. , Yang Z. , Zhu Z. 1995 A new k-ε eddy-viscosity model for high Reynolds number turbulent flows. Computers & Fluids 24 (3) 227 - 238
Silva L.F.L.R. , Lage P.L.C. 2011 Development and implementation of a polydispersed multiphase flow model in OpenFOAM. Computers and Chemical Engineering 35 (12) 2653 - 2666
van Leer B. 1979 Towards the ultimate conservative difference scheme. Journal of Computational Physics 32 (1) 101 - 136