Below are the first 10 and last 10 pages of uncorrected machine-read text (when available) of this chapter, followed by the top 30 algorithmically extracted key phrases from the chapter as a whole.

Intended to provide our own search engines and external engines with highly rich, chapter-representative searchable text on the opening pages of each chapter.
Because it is UNCORRECTED material, please consider the following text as a useful but insufficient proxy for the authoritative book pages.

Do not use for reproduction, copying, pasting, or reading; exclusively for search engines.

OCR for page 633

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
SESSION 12
LIFTING-SURFACE FLOW: STEADY VISCOUS METHODS

OCR for page 633

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
This page in the original is blank.

OCR for page 633

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
Numerical Simulation of Turbulent Flows Around Hydrofoil
C.C.S.Song and C.Chen
(University of Minnesota, USA)
Abstract
The numerical solution to the unsteady compressible flows at small Mach number has always run into trouble due to the disparity between the sound speed and convective speed. The equation of unsteady incompressible flows are even more difficult to solve numerically than that of compressible flows. More importantly, the incompressible flow equations contain no physical information related to pressure waves which is an essential mechanism of rapidly accelerating flow such as hydraulic transients and acoustics. The compressible hydrodynamic equations have been developed to remedy the mathematical and physical shortcomings of the incompressible flow approach. The large eddy simulation approach based on the compressible hydrodynamic equations has been applied to flow around a two dimensional hydrofoil and tip vortex flow around a three dimensional hydrofoil. The simulated mean flow, as well as the turbulent correlations in the two dimensional case have been compared with experimental data. The tip vortex trajectory and the structure of tip vortex core have been studied.
Nomenclature
A+
constant
As
projected surface area
a0
sound speed
Cl
lift coefficient
C0
base chord length
Cs
Smagorinsky constant
D
damping coefficient
L
lift
M
Mach number
n
normal coordinate of wall
p
pressure
p0
reference pressure
Rc
radius of curvature
Re
Reynolds number
Rex
Reynolds number based on x
Rθ
Reynolds number based on momentum thickness
r
radius in polar coordinate
rc
core radius
Sij
strain rate tensor
St
Strouhal number
t
time
U
reference velocity
u
resolved velocity components
ue
velocity at outer flow
unresolved velocity components
ut
tangential velocity
shear velocity
Vt
mean tangential velocity
x,y,z
coordinates
yw
distance from the wall
Δ, Δi
length scales
Г
circulation
ν
molecular viscosity
νt
turbulent viscosity
ω
mean vorticity
ωm
local maximum mean vorticity
ρ
density
ρ0
reference density
θ
momentum thickness
1. INTRODUCTION
The compressible hydrodynamics concept has been commonly used to treat hydrodynamic transient. However, in most cases, only one-dimensional situations have been considered. For three dimensional flow phenomena, the commonly used equations are either compressible Navier-Stokes equation for relatively high Mach number flows, or incompressible Navier-Stokes equations for hydrodynamic flows. The incompressible hydrodynamic equations have been

OCR for page 633

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
widely used in treating fluid mechanics problems with water as the flow medium. But the incompressible assumption changes the nature of the governing equations from that of hyperbolic type to that of elliptic type. Hence, the traditional numerical method for solving hyperbolic equations cannot be applied. Many difficulties have been encountered in solving the incompressible equations.
In order to obtain the hyperbolic property of the incompressible equation, Chorin [1] has proposed an artificial compressibility method. With this method, one can obtain the steady state solution of incompressible flow in the same way as one solves the compressible Navier-Stokes equations. However, for unsteady hydrodynamic problems, this method does not capture the essential effect of compressibility on fluctuating pressure or flow noise. In most unsteady flows, the weak compressibility of the fluid can contribute a great deal to the pressure field, hence the whole flow phenomena.
The three dimensional compressible hydrodynamic equations, which previously were called the weakly compressible model, has been established by Song and Yuan [2]. This set of equations not only retains the hyperbolic property, but also preserves the hydrodynamic compressibility effect in the equations. By using the property of the small change in density in the case of low Mach number flows, the dynamic part of the flow field is decoupled from the energy equation. But the compressibility effect on the pressure field still remains. Compressible hydrodynamic equations have been successfully used in solving various unsteady, low Mach number flow problems [3, 4].
Flow around foils is a class of common and important phenomena in fluid mechanics. Its significance in engineering practice is evident. Tip vortex as one of the flow pattern in three dimensional cases has not been understood fully; however, its effects on the foil performance is essential. It is theoretically interesting and practically important to understand the principle, origin, and development process of tip vortex.
It is extremely difficult to fully resolve large Reynolds number flows because they contain a very large range of different size eddies. The large eddies simulation (LES) approach, by which large eddies are resolved and small eddies are modeled, has been rapidly developed. The current paper presents numerical simulations of turbulent flows around a two dimensional and a three dimensional hydrofoil by solving the compressible hydrodynamic equations with the large eddy simulation method. The tip vortex flow phenomena will be analyzed as well.
2. GOVERNING EQUATIONS FOR LES
2.1 Compressible Hydrodynamic Equations
To implement the large eddy simulation, one needs to apply the cell volume integration over the governing equations, namely, the conservation of mass and momentum equations. A simple volume average method, by which the resolved quantity is treated as a constant over each finite volume, has been used. This is similar to the method introduced by Deardorff [5] and later extended by Schumann [6]. The volume averaged conservation equations have the similar form as the Reynolds average equations. One extra term is generated from the nonlinear term of the momentum equation. It has been referred as the sub-grid scale Reynolds stress which is analogous to the time averaged Reynolds stress.
Customarily, the sub-grid scale Reynolds stresses are decomposed into the sum of a tracefree tensor and a diagonal tensor. The diagonal tensor part is brought out in a similar way as the dynamic pressure p from the viscosity stress tensor. It is usually combined with pressure as a modified pressure, because they have similar properties, except one is from the molecular viscosity stress tensor, the other is from the unresolved turbulent stress tensor. The modified pressure properly represents the resolved pressure field. The trace-free part of the sub-grid scale Reynolds stress tensor need to be modeled in order to close the set of equations. The model will be discussed in more detail later.
The essential point for deriving compressible hydrodynamic equations is to use the equation of state appropriate to liquids or small Mach number, for which density change in the flow is small. Consider a barotropic flow with low Mach number, one can approximate the equation of state as:
(1)
Where p0 and ρ0 are reference pressure and density. a0 is the sound speed, which can be considered as constant. This expression has been widely used in hydraulic transient analysis.
From the above equation, one can solve for

OCR for page 633

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
ρ and substitute it into the averaged conservation equation of mass and hence get:
(2)
By dimensional analysis, one can show that the second term is of order M2 while the first term is of order M2St. In the case of small Mach number flows, the second term can be ignored. But since the first term is related to the Strouhal number, this term cannot be neglected, especially, when flow is highly time dependent and the Strouhal number is large. It is essential to retain the first term in the continuity equation for flows involving rapid acceleration such as hydraulic transient and hydroacoustics. It also keeps the whole equation set to be of the hyperbolic type. It should also be pointed out that the third term in Eq. (2) is the noise production term.
Similarly, one can substitute the equation of state into the conservation equation of momentum. By considering the density change to be small, one can get the final equation set as following:
(3)
where is stress tensor,
Eq. (2) and (3) are called compressible hydrodynamic equations. There are four unknowns and four equations. Obviously, this set of equation is closed except that the sub-grid scale turbulent stress need to be modeled. The energy equation has been decoupled. And the equation set has the hyperbolic property and it retains the compressibility effect to the flow field through the continuity equation, Eq.(2).
2.2 SGS Model
The Smagorinsky model is a simple and widely used model for modeling the unresolved sub-grid scale Reynolds stresses. In analogy to the definition for viscosity stresses, the sub-grid scale Reynolds stresses are assumed to be proportional to strain rate of the resolved flow field [7], i.e.
(4)
where νt is the unresolved turbulent viscosity, or the sub-grid scale diffusivity. νt depends on the length scale of the cell and the strain of the resolved flow field in the following way:
(5)
Where Δ represents the filter width, which is taken to be the length scale of the finite volume. There are several ways to decide Δ. According to Bardina et al. [8], a better choice might be:
whereΔi is the three dimensions of the finite volume.
Cs is the Smagorinsky constant. Many different values for different flow situations have been used. For homogeneous isotropic turbulence, Lilly [9] determined that Cs≈0.23. Deardorff [5] used Cs=0.1 for his simulation of turbulent channel flow. And several simulations [10, 11] have used values between 0.1 and 0.23. It is supposed this constant is a universal constant and independent of grid size. But this seems not the case. The simulations that have been carried out by the authors suggested that this constant needs to be adjusted somehow for different mesh systems and flow situations. For the two dimensional simulation to be presented herein, better results have been obtained by using Cs=0.14. And for the primary tip vortex simulation, the tip vortex trajectory is insensitive to this constant.
For turbulence near a solid wall, Eq. (5) has been modified in a similar way as proposed by Moin and Kim [12].
(6)
Where D is a non-dimensional damping factor. Van Driet exponential damping function has been used,
D=1−exp(−y+/A+)

OCR for page 633

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
It is easy to see, this damping factor will only modify the turbulent viscosity near the solid wall.
3. COMPUTATIONAL CONSIDERATIONS
The size of the largest eddies and the size of the smallest eddies are extremely far apart in turbulent flows. And their difference increases with increasing Reynolds number of the flow. It is very difficult to solve the equations governing the turbulent flows, especially in the case of large Reynolds number flows. In order to gain maximum accuracy with available computer resources, nonuniform mesh system and multi-mesh system has been used.
As the boundary layer on the foil surface and the wake behind the foil are the primary objects in the simulation of flow around two dimensional foil, the computational domain is deliberately divided into two regions. The fine mesh region encloses the foil and covers the wake behind the foil where the gradients are relatively large. The coarse mesh region covers the rest of the computational domain from the lower channel wall to upper channel wall. For the simulation of flow around three dimensional foil, tip vortex is the major objective at this stage of the research. Hence fine mesh has been arranged to resolve the tip part of the foil and the tip vortex behind it. This multi-mesh system has been very beneficial in a sense of using the computer time and memory efficiently.
The artificial interface condition between the fine mesh region and the coarse mesh region is given special attention. For boundary points of the coarse mesh, the primary variables are directly taken from the the corresponding fine mesh region. For boundary points of the fine mesh region, the velocity values and the pressure gradients have been matched to the coarse mesh region. Matching pressure gradient instead of matching pressure itself has enabled the unsteady pressure waves generated by the foil to propagate outward properly. Otherwise, a pressure wave will reflect back and accumulate large errors in the inner zone and hence cause the computation to be broken down.
Wall boundary conditions are very important for the computation. Special attention should be paid to high Reynolds number flows, since the mesh size usually is not fine enough to resolve the details of the boundary layer. Appropriate partial slip boundary conditions have to be used. Many [3, 6, 13] have used the wall function as a medium to obtain the approximation for the phantom point velocity values. However, the associated assumption of fully developed turbulent boundary layer has limited its usage. Especially it cannot be applied for flow that has laminar boundary layer, transition, as well as turbulent boundary layer. A new way to impose the partial slip condition has been developed for the present simulation in the case of two dimensional flow. By assuming the boundary layer is quasi-steady during a certain time period, say, the time interval of the outer flow calculation increment, one can use the steady boundary layer equations by ignoring the local acceleration term, .
In the current simulation, Crank-Nicolson's fully implicit second order scheme is used to solve the boundary layer equations. The outer flow conditions are obtained by inputting the pressure field near the wall from the outer flow computation. In the leading part of the foil, the laminar boundary layer is assumed. The transition starting point is determined by laminar boundary layer separation or according to the following criteria [14],
where
θ is the momentum thickness of the boundary layer.
Starting from whichever comes first, the turbulent boundary layer equation will be applied. For turbulent boundary layer calculation, Prandtl's mixing length model is used.
After the boundary layer separation, non-slip velocity boundary condition has been used, since after separation, the boundary layer is fairly thick, and the velocity gradient normal to the wall has become much smaller compared with the boundary layer before the separation. Hence, the non-slip boundary condition will be appropriate.
The quasi-steady boundary layer calculation provided the information for imposing the partial slip boundary conditions in the outer flow. It is a novel combination of boundary layer computation and outer flow simulation. For unsteady

OCR for page 633

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
large eddy simulation, it certainly can improve the accuracy in reflecting the boundary layer effect to the outer flow and flow in the wake.
The pressure condition on the solid walls has been specified in the following way,
(7)
Where Rc is the radius of curvature of the wall. This is obvious, since the curved wall will generate a centrifugal force, which will result in a pressure gradient normal to the wall.
Boundary layers along the channel walls are not so interesting for present study; in addition, these boundary layers will not affect the foil wall boundary layer and the wake region. The channel wall has been simply treated with full slip velocity and zero gradient pressure conditions.
The well-known MacCormack predictor-corrector numerical scheme has been used in solving the compressible hydrodynamic equations for large eddy simulation. The scheme is of second order accuracy in time and space.
4. NUMERICAL RESULTS AND DISCUSSIONS
4.1 Flow Around 2-D Foil
A two dimensional foil that has been tested in a water channel by Huang [15] at National Taiwan University has been simulated. The Reynolds number of the flow is 1.6×106. And the angle of attack is −5.9.
Fig. 1 shows two instantaneous vorticity fields near the trailing edge of the foil and in the wake. One may note that at pressure side of the foil, a very thin boundary layer is formed, and the vorticity in this layer has a positive sign. On the other hand a relatively thick boundary layer exists at the suction side of the foil. The vorticity is negative in this layer. Apparently, there is no visible separation for the simulated flow. However, the vorticity with opposite sign from two sides of the foil meet at the trailing edge and form a vortex sheet, which is very unstable. Vortices roll-up and negative and positive vortices alternatively shed downstream. A more vivid picture of the
Fig. 1 Instantaneous Vorticity Contours

OCR for page 633

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
vortex shedding can be seen through animations of color coded pressure and vorticity fields.
The mean flow field around the foil body is shown in Fig. 2. The calculated values are compared with the experimental data. It can be seen they agree very well. The mean flow profiles in the boundary layer are plotted against the experimental measurement in Fig. 3. Again, fairly good agreement has been reached. Mean flow velocity
Fig. 2 Simulated and Measured Mean Flow Ux around the Foil
Fig. 3 Simulated and Measured Mean Boundary Layer Profiles
vectors are shown in Figs. 4 and 5 for simulated values and experiment values respectively. The agreement between simulated values and experimental data is quite good. It can be concluded that the mean velocity field can be reproduced by the numerical simulations very well.
As discussed above, the turbulence has been divided into the resolved large scale part and the unresolved sub-grid scale part. The former is simulated while the latter is modeled through the resolved flow field. Figs.6(a) and (b) provide the
Fig. 4 Simulated Velocity Vectors in the Wake
Fig. 5 Measured Velocity Vectors in the Wake
turbulent intensity for these two parts in several stations in the wake. And the sum of the two parts is plotted in Fig. 6(c). It is interesting to note that in the wake region near the trailing edge, the main contribution is from the modeled part, while in the wake region far from the trailing edge, the main contribution is from the resolved large eddies. This is consistent with the general understanding that eddies are larger in the far field but they are smaller near the trailing edge. The pattern that eddy sizes change from small to large as flow goes downstream can be easily seen on Fig. 1. In Fig. 6(d), the measured turbulent intensity is shown. Comparing Fig.6(c) and Fig. 6(d), one can see that the shape is quite similar, mainly in the far wake region. But in the near wake region, the sharp peak right after the trailing edge as shown in Fig./, 6(d) is lost in the simulation. Evidently a very thin shear layer after

OCR for page 633

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
the trailing edge produces the strong and fine turbulent eddies. Relatively poor agreement in the near field may be due to the three dimensionality of small scale and insufficient grid refinement.
Fig. 6 Comparison of Turbulent Intensity in the Wake
4.2. Tip Vortex Flow
A three dimensional foil with a modified NACA 4215 section and an elliptic planform has been studied. Extensive experimental research for tip vortex cavitation with this foil has been conducted [16, 17]. Tip vortex is the focus in the present numerical simulation. The coordinate axes are so chosen that x is along the main flow direction, y is normal to the plane of the foil, positive pointing from pressure side to suction side, and z is in the span-wise direction from the foil base to the tip of the foil. The effective angle of attack is about 14º. And the Reynolds number of the simulation is 7.0×105. However, since the mesh size are not fine enough and the boundary layer on the foil surface is not resolved well in the current three dimensional simulation.
The lift coefficient is defined as Cl= , where L is the total lift; As is the projected area of the foil. Calculated lift coefficient is 0.72, which is very close to the experimental result, 0.72∼0.73 [ 17].
Fig. 7 shows the velocity vectors projected on a cross-section normal to the flow direction. The swirl around the tip vortex is quite apparent. The interaction between wake and tip vortex can be observed by noting the sudden changing of flow direction below the swirl. The corresponding vorticity distribution (x-component) is

OCR for page 633

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
shown in Fig. 8. The vorticity contour lines form nearly concentric circles indicating the intensity
Fig. 7: Projected Velocity Vectors at x/C0= 1.06 Cross-Section
Fig. 8 Vorticity Contours of x-Component at x/C0=1.06 Cross-Section
increases towards the vortex core. Vorticity contours on a surface parallel to the x-z plane and through the axis of the tip vortex are shown in Fig. 9.
The simulated tip vortex trajectory
Fig. 9 Vorticity Magnitude at y/C0=0.0 Cross-Section
Fig. 10 Tip Vortex Trajectory
is plotted and compared with the measurement in Fig. 10. The agreement is very good. It is worth pointing out that the current simulation has not tried to resolve the boundary layers on the foil surface. It is believed that the boundary layer on the foil surface germinates the tip vortex, but the above result seems to imply that the tip vortex trajectory is independent of the boundary layer. This result agrees with the measurement [16] that the tip vortex trajectory is insensitive to the Reynolds number and the angle of attack.
The pressure coefficient in the tip vortex core as a function of x is plotted in Fig. 11. It can be seen that the pressure reaches the lowest point at the tip where the tip vortex originates. The pressure increases in the downstream direction. This is understandable because the tip vortex core is very small at the beginning and grows in the downstream direction due to diffusion. Increasing core size means weakening vorticity and increasing pressure.

OCR for page 633

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
Fig. 11 Minimum Pressure Coefficient at Tip Vortex Core
It is a difficult task to measure the full details of the structure of tip vortex. On the other hand theoretical analysis of tip vortex has always required an assumption that the tip vortex is a Rankine vortex. The Rankine vortex model ignores two aspects of the real fluid flow phenomena. Firstly, according to the Rankine vortex model, the vorticity field should be axis symmetric about the vortex core. However, the wake behind the foil always disturbs this symmetry. The region affected by the wake becomes larger as the wake spreads wider along the streamwise direction. Measurements engaging the assumption of symmetric vorticity field could lead to errors in their data since the tangential velocity profiles around the tip vortex core are quite different depending on how the wake affects them. Secondly, the Rankine vortex model assumes that, the flow field around tip vortex is divided into rotational core region and irrotational outer region. And the vorticity is concentrated inside the core region. This cannot be true for the real flow situation because of the diffusivity of the fluid and turbulence. Again the wake plays an important role in distributing the vorticity in the flow field.
Numerical diffusion can distort the results of numerical simulations. But numerical simulations can provide at least a qualitative structure of the real tip vortex. Since the tip vortex is not symmetric and the vorticity is not just concentrated inside the core region, the traditional definition for core radius needs to be reconsidered. In the present study, the circulation Г about the core is calculated first as following,
Fig. 12 Nondimensional Circulation around Tip Vortex Core
The variations of Г as function of radius from the core at several stations along streamwise direction have been plotted in Fig. 12. Note that Г is nondimensionalized by taking theoretical value of as the reference value so that, in the ideal case, all curves in Fig. 12 would approach 1 asymptotically. The asymptotic values for the last four curves is about 1.05. This slight difference from 1 could be caused by the fact that the flow simulated here is bounded by the walls. By comparing the curves for different stations, one could see that the asymptotic value starts with about 0.75 at the tip and gradually increases to the maximum value of 1.05 at about x/C0=0.5. This means that the vorticity in the z-direction is gradually changing the direction into the vorticity in the x-direction, i.e., the tip vortex. In other word, the circulation contributed to the lift is gradually converting into the circulation of tip vortex. This process finishes at about x/C0=0.5 position. The steep gradient near the core shows that the vorticity is very strong near the core, and becomes weak as it goes far from the core.
The mean tangential velocity can be calculated based on the circulation calculated above. The mean tangential velocity at one of the x-station has been shown in Fig. 13. The corresponding tangential velocity profile of the Rankine vortex is compared in the same figure. As can be seen, in the core region, the two are quite close. But outside the core, they are very different. No-

OCR for page 633

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
In the near wake station, , the comparison between the cross-flow velocities of the two solutions is similar to the one described for the wing station. The differences of the two solutions close to the external boundary of the reference case are larger than in the wing station, but the effect in the tip vortex region is still very small.
Figure 9: Pressure distribution at the trailing edge computed from the upper and lower surfaces of the wing for solutions obtained with different locations of the external boundary.
The pressure discontinuity at the trailing edge is illustrated in figure 9 for the calculations with the two different locations of the external boundary. The larger pressure discontinuities are obtained in the tip vortex region, where a small pressure oscillation is obtained on the lower surface for both solutions. A small difference between the pressure level of the two solutions is obtained, but the pressure discontinuity is the same in the two calculations. These results indicate that the pressure discontinuity is mainly determined by the difference in circulation of the outer flow in the potential and viscous flow calculations, which is not affected by the change in the location of the external boundary.
The comparison of the results of the calculations performed with different locations of the external boundary showed that the influence of the approximate boundary conditions imposed at the external boundary on the solution in the tip vortex region are small. These results suggest that the present zonal approach is acceptable for tip vortex calculations, because it reduces substantially the size of the region to be discretized, compared with an alternative external boundary at a location where undisturbed flow conditions could be imposed. Further investigations will be required to address the influence of the viscous/inviscid interaction.
4.3 Specification of the initial velocity profiles
At the inlet station two types of boundary conditions were tested : prescribed velocity profiles and prescribed streamwise velocity gradients. It should be mentioned the the Neumann condition was applied as a Dirichlet condition, updated in the course of the solution process. So in the case of prescribed streamwise velocity gradient the velocity profiles at the inlet boundary are obtained iteratively. This makes the convergence usually slower than in the case of prescribed velocity profiles. However, the specification of velocity profiles at the inlet station may be difficult and a poor approximation of the initial velocity profiles can cause divergence of the solution close to the inlet station.
Our approach to the specification of the initial velocity profiles is based on simple assumptions. The velocity profiles are derived from standard boundary layer profiles, [22], defined by the local momentum thickness, θ, and skin friction coefficient, Cf. The values of θ and Cf at the inner boundary are obtained from a 2-D calculation, [27], and the spanwise variation of θ and Cf is obtained by a Hermite interpolation of the predicted 2-D values on the upper and lower sur-

OCR for page 633

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
faces at the inner boundary. In the potential flow region the velocity profiles are obtained from a potential flow solution, [21]. Due to numerical difficulties with the source-based panel method, it is not possible to calculate the potential flow solution at very small distances to the wing surface and so we have calculated the velocity components along a line η=constant at a sufficient distance from the wing surface. In the region between this line and the viscous region thickness, determined by the specified θ and Cf, we assume that the potential flow velocity components do not change in the η direction.
The streamwise velocity gradient is specified applying the procedure described above to the first two stations of the grid and using a first order approximation to the derivatives. With this weaker boundary condition of the Neumann type, the velocity profiles at the inlet boundary are allowed to adjust to the local flow characteristics. As mentioned before, the streamwise velocity gradient is not specified implicitly in the calculation. The velocity profiles at the inlet boundary are updated in the course of the iteration process to obtain the specified streamwise velocity gradient. This makes the convergence rate of the solution lower and highly time consuming because the initial profiles change between consecutive sweeps. In the present solution the inlet velocity profiles were updated only when the maximum pressure difference between consecutive sweeps at the second streamwise station was ΔCp>2.0×10−2.
The flow of the reference case was calculated with both types of inlet conditions. The results obtained with specified velocity profiles at the inlet boundary exhibit a rapid change of the pressure field in the initial streamwise stations. A clear visualization of the disturbances caused by the specification of the inlet velocity profiles is presented in figure 10, where the pressure field at the first calculated station, , of the two solutions is compared. The pressure field obtained with specified inlet velocity profiles is clearly disturbed by the specification of the initial velocity profiles. With specified streamwise velocity gradient a smooth pressure field is obtained at .
We found that the disturbances in the solution obtained with specified velocity profiles are restricted to the three initial stations and so it is important to know the effect of the inlet boundary conditions in the generation and development of the tip vortex. The calculated pressure fields of the two solutions are compared in figure 11,
Figure 10: Comparison of the pressure fields at the first streamwise station for solutions obtained with different inlet boundary conditions.
for the streamwise stations at and .
In the wing station, , the minimum pressure coefficient is 0.05 lower in the solution obtained with Neumann boundary conditions at the inlet station. The differences are smaller than the ones obtained at the inlet stations, but the effect of the different inlet boundary conditions is still visible in this region. However, the vortex location is virtually the same in the two calculations. In the wake stations, in figure

OCR for page 633

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
Figure 11: Comparison of the pressure fields at two different streamwise stations for solutions obtained with different inlet boundary conditions.
11, the comparison of the two solutions is similar to the one described above for the wing station. The differences between the two solutions have become smaller than in the wing station; the tip vortex location is again hardly affected.
4.3.1 Effect of the Tip Geometry
It is interesting to know whether the behaviour of the solution process for the two types of inlet boundary conditions is the same for other tip geometries. Therefore the flow at the tip of a rectangular wing with a squared tip was also calculated. The wing aspect ratio (4) and airfoil section (NACA 0015) are the same as in the reference case. The same angle of attack and Reynolds number were used, respectively 6.5 degrees and 8.5×105. In the previous test cases the tip geometry is smooth. In a wing with a squared tip the geometry has two discontinuities in the tangential derivatives. The local features of the flow close to the geometry discontinuities make the specification of the inlet velocity profiles more difficult in this case. In the present specification of inlet velocity profiles no special treatment was included to account for the geometry discontinuities, the simple procedure described before was applied without any modifications and so the inlet conditions will locally be highly inaccurate. The grid used in these calculations is similar to the grid of the reference case illustrated in figure 3. The location of the boundaries and the streamwise stepsize are equal to the ones used in the reference case.
In the calculation performed with specified velocity profiles at the inlet boundary the pressure field did not converge in the initial stations of the grid. The pressure field showed a non-smooth behaviour in these initial stations and several jumps in (ΔCp)max are obtained in the first streamwise station downstream of the inlet plane. Nevertheless in the last sweeps performed, the pressure differences between consecutive sweeps downstream of the initial streamwise stations are three orders of magnitude smaller than the pressure differences close to the inlet boundary.
The same flow was calculated with specified streamwise velocity gradient at the inlet boundary. In this solution the inlet velocity profiles were updated only when the maximum pressure difference between consecutive sweeps at the second streamwise station was ΔCp>2.0×10−2. With this type of boundary conditions it is possible to obtain a converged pressure field in the whole computation domain. The imposition of a weaker boundary condition at the inlet boundary removed the jumps of (ΔCp)max at the inlet boundary obtained with specified velocity profiles.
The pressure field at the first calculated sta-

OCR for page 633

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
Figure 12: Comparison of the pressure fields at the first streamwise station for solutions obtained with different inlet boundary conditions.
tion, , of the two solutions is compared in figure 12. The disturbances in the pressure field close to the inlet boundary are evident and a particularly non-smooth solution is obtained close to the geometry discontinuities. A smoother solution is obtained at with specified streamwise velocity gradient.
The effect of the different inlet boundary conditions on the tip vortex generation and development is similar to the one obtained on a wing with a rounded tip.
These results show that the specification of inaccurate velocity profiles at the inlet boundary can deteriorate the convergence of the method close to the inlet boundary. The results obtained with streamwise velocity gradient at the inlet boundaries did not reveal any convergence difficulties of the method and the use of weaker boundary conditions enabled the adjustment of the solution to the local flow characteristics. This suggests that innaccuracies in the specification of the streamwise velocity gradient are less damaging to the solution process than inaccuracies in the specification of inlet velocity profiles. However, with Neumann boundary conditions at the inlet boundary the method becomes more time consuming than with specified velocity profiles.
4.4 Comparison with Experimental Results
The results obtained with the present method are compared with experimental results available in the literature, and with the results of the method of Govidan et al., [2], to evaluate its performance. The test cases include two wings of rectangular planform with a squared and a rounded tip.
4.4.1 Wing with Rounded Tip
The first test case is a rectangular wing of aspect ratio 6 with a NACA 0012 airfoil section and a rounded tip. The purpose of this test case is to compare the results of the present method with the results of the method of Govidan et al., [2]. The latter also uses a space-marching process to obtain the solution, but with a single sweep through the computational domain in the streamwise direction, because a full parabolisation of the Navier-Stokes equations is used, the streamwise pressure gradient is obtained from an inviscid flow solution. Although a different solution method is adopted by Govidan et al., the physical approximations assumed in the governing equations are the same as the ones used in the present method with the exception of the streamwise pressure gradient. The comparison of the results of the two methods enables the evaluation of the streamwise pressure gradient influence on tip vortex calculations. The calculations were performed for the conditions given in [2] : Angle of attack of 6.18 degrees and Reynolds number based on the chord length of 7.52×105.

OCR for page 633

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
Figure 13: Comparison between the predicted pressure distribution on the wing surface and experimental results for the flow at the tip of a rectangular wing, with an angle of attack of 6.18 degrees and a Reynolds number of 7.52×105.
The location of the boundaries of the computation domain were identical to the ones used by Govidan et al. in [2], with the exception of the outlet boundary. The results presented in [2] include only the wing region, the outlet station is 6% of the chord upstream of the trailing edge. In the method of Govidan et al., the flow solution at a transverse plane depends only on upstream occurrences and so there are no streamwise boundary conditions required at the outlet plane. In the present method, the elliptic character of the equations is retained in the pressure field and so a downstream pressure boundary condition is required at the outlet station. Close to the trailing edge of the wing, it is not possible to impose a physically acceptable boundary condition and so the outlet station has to be moved into the far wake. In the present grid, the outlet station is located three chords downstream of the trailing edge. The inlet boundary is 15% of the chord downstream of the leading edge. Like in the calculations of Govidan et al., the inlet boundary conditions of the present calculations consist of specified velocity profiles. The inner boundary, inlet and outlet of the cross-flow, is placed at a distance of 40% of the chord length inboard of the tip of the wing. The distance between the external boundary and the wing is 30% of the chord length at the station of maximum thickness of the wing.
A number of 161 streamwise stations was used, including 81 on the wing and 80 on the wake. A 40×47 2-D grid was generated on each streamwise station11.
In figure 13 the calculated pressure distribution on the wing surface is compared with the predictions of Govidan et al., [2], and the experimental results of Gray et al. [28], given in [2]. The present results show better agreement with the experimental results than the predictions of Govidan et al. Notice the mismatch of predicted and experimental isobars close to the leading edge in the results of Govidan et al. A significant improvement is obtained in the prediction of the pressure distribution close to the tip. We note that Govidan et al., [2], used a 120×40×47 grid to discretize the viscous region between 15% and 94% of the chord, while the present results were obtained with an equal number of nodes in the transverse sections but with only 81 stations between 15% of the chord and the trailing edge.
4.5 Wing with Squared Tip
Experimental results including mean velocity LDV measurements have been reported recently by Falcão de Campos et al., [29], for a rectangular wing of aspect ratio 4 with a NACA 0015 airfoil section and a squared tip. The experimental investigations included flows at three different angles of attack, 4.5, 6.5 and 8.5 degrees at a Reynolds number of 8.5×105. These flows have been numerically simulated with the present method.
The boundaries of the viscous flow region were placed in the same locations as in the reference case of section 4. The grid has 145 streamwise stations, including 81 on the wing and 64 on the wake. In the transverse sections of each streamwise station a 51×47 2-D grid was generated. The calculations were performed with Neumann boundary conditions at the inlet station, i.e. the streamwise velocity gradient is specified.
The experimental measurements were performed on planes perpendicular to the undisturbed flow and the measuring plane does not coincide with a grid plane. To compare the present
11
The grid is similar to the one used in the reference test case illustrated in figure 3.

OCR for page 633

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
Figure 14: Comparison between the present predictions and the experimental results for the flow at the tip of a rectangular wing with a squared tip. Angle of attack 6.5 degrees and Reynolds number of 8.5×105.
results with the experimental data the velocity components at the experimental measuring plane are obtained from the values at the grid nodes using linear interpolation along the ξ lines. Furthermore, the comparisons are performed in the cartesian coordinate system of the experimental investigations, which has the x* direction aligned with the undisturbed flow, U∞. In the present comparisons the origin of the coordinate system (x*,y*,z*), is placed at the tip of the trailing edge, (x*=0, y*=0, z*=0) ≡ (x=1.0,y= 0.0,z=2.0).

OCR for page 633

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
The numerical predictions of the transverse velocity fields at and are compared with the experimental results in figure 14, for the flow with an angle of attack of 6.5 degrees. An excellent correlation of the tip vortex location is obtained in both planes, but in both cases the maximum transverse velocities are underpredicted in the calculations. Upon inspection it may be seen that a more diffused vortex is obtained in the calculations. This is particularly evident at the plane , where a nearly axisymmetric flowfield is obtained in the vicinity of the vortex core. However, the asymmetry of the vortex close to the trailing edge is well predicted and the flowfield away from the vortex center shows a very good correlation with the experimental results. The location of the vortical wake in the station close to the trailing edge, , is also well predicted. Furthermore, although the maximum values of the tranverse velocity components are larger in the experimental results, the downstream decrease of in the predictions is very close to the experimental one.
The turbulence model used in this method does not account for streamline curvature effects. It is well documented in the literature that in the vortex core the turbulent diffusion should be negligible, see for example [30]. In the present turbulence model the eddy viscosity is related to the displacement thickness, which is determined along the grid η lines. The streamwise velocity defect in the tip vortex region causes an increase of the local displacement thickness and consequently of the eddy viscosity. This means that with the present model the eddy viscosity exhibits a maximum at the tip vortex region, which is in total contradiction with the physical characteristics of the flow. It is thought that the large values of the eddy viscosity in the vortex core may be reponsible for the overprediction of the vortex core size in the present calculations.
The overprediction of the vortex core size has a large effect on the calculated pressure in the vortex core. The minimum pressure coefficient in the vortex core is approximately proportional to the square of the inverse of the core size, which means that an accurate prediction of the core size is essential for a good prediction of the pressure field.
The comparison between predictions and experimental results for the flows with angles of attack of 4.5 and 8.5 degrees is similar to the one obtained with an angle of attack of 6.5 degrees, [12]. The correlation between the vortex location of the predictions and experimental results is still excellent, but the vortex core size is again overpredicted. The location of the vortical wake close to the trailing edge is in good agreement with the experimental results for both angles of attack and the flowfield away from the vortex center is in excellent agreement with the experimental results.
5 Conclusions
A numerical method for the calculation of the steady, incompressible, viscous flow at the tip of wings was developed. The effect of the choice of the extent of the computation domain on the solution of tip vortex flows has been investigated by changing the location of the boundaries. The results showed that the present zonal approach to the calculation of tip vortex flows is efficient, because a large part of the flow field does not need to be discretized and a finer discretization can be obtained in the tip vortex region for a given number of grid nodes. However, the imposition of potential flow conditions at the external boundary implies that the viscous/inviscid interaction has to be incorporated in the solution procedure to account for the difference in the external flow circulation associated with the viscous and inviscid solutions respectively. In the present calculations the viscous/inviscid interaction is neglected and so a pressure discontinuity is obtained at the trailing edge.
The specification of appropriate inlet boundary conditions may be difficult. When initial velocity profiles are specified convergence problems in the initial stations may result if a poor estimate of the initial velocity profiles is used. An improved solution can be obtained by specifying the streamwise velocity gradient at the inlet boundary. However, this type of boundary condition reduces the convergence rate of the method, making it more time consuming. The comparison between the solutions obtained with the two types of boundary conditions suggests that the influence of small inaccuracies in the inlet velocity profiles is restricted to a small region close to the inlet boundary.
The present calculations showed that the streamwise pressure gradient of the inviscid flow solution is not a good approximation to the streamwise pressure gradient of the viscous flow in the tip vortex region. This was revealed in the comparison with experimental results, [28], of the present method and the method of Govidan et al.,

OCR for page 633

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
[2], which imposes the streamwise pressure gradient from inviscid flow calculations.
A good correlation between numerical predictions and experimental results was obtained for the calculation of the flow at the tip of a rectangular wing with a squared tip. However, the size of the tip vortex core is overpredicted in the present calculations. The turbulence model is likely to be responsible for the excess of diffusion in the tip vortex core, because the physical characteristics of the local flow are not well reproduced in the present model. The overprediction of the vortex core size still precludes an accurate prediction of the suction peak in the vortex core.
Acknowledgement
The work was carried out while the first author was on leave from Institute Superior Técnico and supported by grants of the Comissão Permanente INVOTAN of the Junta Nacional de Investigação Cientifica and of the SCIENCE PLAN of the EUROPEAN COMMUNITIES which are gratefully acknowledged. The support of the Stichting Nationale Computer Faciliteiten in the present calculations is also gratefully acknowledged.
References
[1] Govidan T.R., Levy R., Shamroth S.J. — Computation of The Tip Vortex Generation Process for Ship Propeller Blades—4th International Conference on Ship Hydrodynamics, Washington D.C., 1985.
[2] Govidan T.R., Jong F.J., Levy R., Shamroth S.J. —Validation of a Forward Marching Procedure to Compute The Tip Vortex Generation for Ship Propeller Blades— 17th Symposium on Naval Hydrodynamics, August-September 1988, The Hague, The Netherlands.
[3] Briley W.R., McDonald H. —Three-Dimensional Viscous Flow with Large Secondary Velocity. —Journal of Fluid Mechanics, Vol. 144, 1984, pp. 47–77.
[4] Rubin S.G., Tannehill J.C. —Parabolized/Reduced Navier-Stokes Computational Techniques. —Annual Review of Fluid Mechanics, Vol. 24, 1992, pp. 117– 144.
[5] Rubin S.G. —Incompressible Navier-Stokes and Parabolised Navier-Stokes Formulations and Computational Techniques . — Computational Methods in Viscous Flows, Vol. 3 in the series Recent Advances in Numerical Methods in Fluids (ed. Habashi) Pineridge Press, 1984.
[6] Rubin S.G. —Global Relaxation Procedure for a Reduced Form of the Navier-Stokes Equations. —Proceedings of the 9th—International Conference on Numerical Methods in Fluid Dynamics, Lecture Notes in Physics, Vol. 218, Springer-Verlag, 1985, pp. 62–71.
[7] Raven H.C., Hoekstra M. —A Parabolised Navier-Stokes Solution Method for Ship Stern Flow Calculations. —2th International Symposium on Ship Viscous Resistance, Goteborg Sweden, March 1985.
[8] Hoekstra M., Raven H.C. —Application of a Parabolised Navier-Stokes Solution System to Ship Stern Flow Computation . —Osaka International Colloquium on Ship Viscous Flow, Osaka Japan, October 1985.
[9] Hoekstra M., Raven H.C. —Ship Boundary Layer and Wake Calculation with a Parabolised Navier-Stokes Solution System. —4th International Conference on Numerical Ship Hydrodynamics”, Washington D.C., 1985.
[10] Hoekstra M. -Recent Developments in a Ship Stern Flow Prediction Code. —5th International Conference on Numerical Ship Hydrodynamics, Hiroshima, September 1989.
[11] Rosenfeld M., Israeli M., Wolfshtein M. — A Method for Solving Three-Dimensional Viscous Incompressible Flows over Slender Bodies—Journal of Computational Physics, Vol. 88, 1990, pp. 255–283.
[12] Eça L. —Numerical Solution of the Parabolised Navier-Stokes Equations for Incompressible Tip Vortex Flows —PhD Thesis, Institute Superior Técnico, Lisbon, March 1993.
[13] Vinokur M. —Conservation Equations of Gasdynamics in Curvilinear Coordinate Systems. —Journal of Computational Physics, Vol. 14 , 1974, pp 105–125.

OCR for page 633

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
[14] Warsi Z.U.A. —Conservation Form of the Navier-Stokes Equations in General Nonsteady Coordinates. —AIAA Journal, Vol. 19, February 1981, pp 240–242.
[15] Zhang H., Camarero R., Kahawita R. — Conservation Form of the Equations of Fluid Dynamics in General Nonsteady Coordinates. —AIAA Journal, Vol. 23, November 1985, pp 1819–1820.
[16] Hoekstra M. —Some Fundamental Aspects of the Computation of Incompressible Flows. —Second Osaka International Colloquium on Ship Viscous Flow, September 1991. Osaka Japan.
[17] Eça L., Hoekstra M. —Discretization of the Parabolised Navier-Stokes Equations. — First European Computational Fluid Dynamics Conference, Brussels, September 1992.
[18] Mynett A.E., Wesseling P., Segal A., Kassels C.G.M. —The ISNaS Incompressible Navier-Stokes Solver: Invariant Discretization. —Applied Scientific Research, Vol. 48, 1991, pp. 175–191.
[19] Segal A., Wesseling P., Van Kan J., Oosterlee C.W., Kassels K. —Invariant Discretization of the Incompressible Navier-Stokes Equations in Boundary-Fitted Coordinates. —International Journal of Numerical Methods in Fluids, Vol. 15, 1992, pp. 411–426.
[20] Sokolnikoff I.S. —Tensor Analysis. —John Wiley & Sons Inc., 1951.
[21] Raven H.C. —Berekening van de potentiaalstroming rond draagvlakken met het programma DAWSON, (in Dutch) ”Calculation of Potential Flow on Lifting Surfaces with the Program Dawson, (Engl. Transl.) —MARIN Report No 50501–1-RD, May 1985.
[22] Hoekstra M. —Generation of Initial Velocity Profiles for Boundary Layer Calculations. — Marin Report No 50028–1-SR, March 1980.
[23] Launder B. —On the Modelling of Turbulent Industrial Flows. —First European Computational Fluid Dynamics Conference, Brussels, September 1992.
[24] Rodi W., Scheuerer G. —Scrutinizing the k–∈ Model Under Adverse Pressure Gradient Conditions. —Journal of Fluids Engineering, Transactions of the ASME, Vol. 108 , June 1986, pp. 174–179.
[25] Cebeci T., Smith A.M.O. —Analysis of Turbulent Boundary Layers. —Academic Press, November 1984.
[26] Cebeci T., Clark R.W., Chang K.C., Halsey N.D., Lee K. —Airfoils with Separations and the Resulting Wakes. —Journal of Fluid Mechanics, Vol. 163, 1986, pp. 323–347.
[27] Eça L.R.C., Falcão de Campos J.A.C. — Analysis of Two-Dimensional Foils Using a Viscous-Inviscid Interaction Method. — International Shipbuilding Progress, (accepted for publication).
[28] Gray R.B., McMahon H.M., Shenoy K.R., Hammer M.L. —Surface Pressure Measurements at Two Tips of a Model Helicopter Rotor in Hover—NASA CR-3281, 1980.
[29] Falcão de Campos J.A.C., George M.F., Mackay M. —Velocity Measurements of the Tip Vortex Flow in the Near- Wake of Hydrofoils. — 1993, Submitted for publication to the Journal of Fluids Engineering.
[30] Bradshaw P. —Effects of Streamline Curvature on Turbulent Flow—AGARDograph 169, AD—768316, 1973.

OCR for page 633

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
DISCUSSION
by Dr. Spyros A.Kinnas
MIT
The authors presented a very interesting approach in solving the tip vortex simulation problem. I would like to ask two questions:
They match a potential flow solution in the inner part of the span to a viscous representation of the outer part. How sensitive are the results on the cut-off location?
In the case of a rectangular planform wing with round geometry in the spanwise direction, were they able to capture a similar vortex roll-up structure along the chord to that of a body of revolution at an angle?
Author's Reply
Our reply to the questions asked by Dr. Kinnas is:
In our present approach, see section 2.1, the boundary conditions in the inner part of the span , imply that cross-stream derivatives, , of the cartesian velocity components and of the pressure are set equal to zero. A potential flow solution is used to specify boundary conditions at the external boundary, without viscous-inviscid interaction. This means that approximate boundary conditions are used on both boundaries. The comparisons between results obtained with different locations of the inner and external boundaries presented in sections 4.2.3 and 4.2.4 show a small influence of the boundaries location on the solution in the tip vortex region.
In this study we did not make any visualization of the limiting streamlines and so we have not yet performed a detailed analysis of the separation pattern at the tip of the wing. The vortex roll-up structure in the case of a round-nosed body of revolution at an angle of attack is strongly dependent on the angle of attack and can be rather complex, as shown for instance in Tobak and Peak, [1]. However, from the analysis of the cross-flow velocities close to the tip, [2], there may be some similarity between the rounded tip wing and the body of revolution at an intermediate angle of attack, where a single dominant vortical structure exists. Such structure has been captured with a similar method in the calculations of Rosenfeld et al., [3], for a prolate spheroid at 10º of incidence.
[1] Tobak, M. and Peake, D.J. Topology of Three-Dimensional Separated Flows— Annual Review of Fluid Mechanics, 1982, Vol. 15, pp. 61–85.
[2] Eca, L. Numerical Solution of the Parabolized Navier-Stokes Equations for Incompressible Tip Vortex Flows, PhD Thesis, Instituto Superior Técnico, Lisbon, March 1993.
[3] Rosenfeld, M., Israeli, M. and Wolfshtein, M. A Method for Solving Three-Dimensional Viscous Incompressible Flows Over Slender Bodies, Journal of Computational Physics, 1990, Vol. 88, pp. 255–283.

OCR for page 633

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
This page in the original is blank.