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 616
lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
AN UNSTEADY THREE-DIMENSIONAL EULER SOLVER COUPLED WITH A CAVITATING PROPELLER ANALYSIS METHOD 616
An Unsteady Three-Dimensional Euler Solver Coupled with a
Cavitating Propeller Analysis Method
J.-K.Choi, S.Kinnas (The University of Texas at Austin, U.S.A.)
ABSTRACT
A fully three-dimensional unsteady Euler solver, based on a finite volume scheme and the pressure correction
method, is developed and applied to the prediction of the effective wake for propellers subject to non-axisymmetric
inflows. The propeller is modeled via unsteady body force terms in the Euler equations. First, the Euler solver is validated
against the analytical solution of actuator disk theory, and the solution in the case of uniform inflow. Then, the unsteady
effective wake is predicted for a three-cycle axial inflow wake, and it is found that the unsteadiness in the effective wake
is small in this case. Finally, the predicted cavity shape and the unsteady velocity field are compared with the
experimental data measured in cavitating propeller experiments which have been performed at the MIT water tunnel.
1 INTRODUCTION
The effective wake is a very crucial issue during the design and assessment of propeller performance, especially in
the presence of blade cavitation, because it is well known that the predicted cavity extent and volume, as well as the
magnitude of the predicted pressure pulses depend strongly on the estimated effective wake inflow.
There has been experimental (Huang et al 1976, Huang & Cox 1977) and theoretical work (Huang & Groves 1980,
Shih 1988) for the prediction of the effective wake by employing steady axisymmetric Euler equations. Later, effective
wake prediction methods using Reynolds Averaged Navier-Stokes (RANS) equations were developed (Stern et al 1988b,
Stern et al 1988a, Kerwin et al 1994, Kerwin et al 1997) for axisymmetric flow applications and (Stern et al 1994) for non-
axisymmetric applications. In the above methods, the propeller is represented by the body force term in RANS equations.
The method using the body force to represent the blades was also used in turbomachinery applications (Damle et al 1997),
where the body force was determined by iteratively adjusting it until the velocity tangency condition on the rotor and
stator blades was satisfied. A different approach which does not need the propeller potential solution is to solve the RANS
equations directly with the no-slip boundary condition on the propeller blade surface. (Chen et al 1994, Stanier 1998)
In a first effort toward the better prediction of the effective wake harmonics, the authors have already developed a
steady three-dimensional Euler solver and successfully applied it to predict the steady non-axisymmetric effective wake
for given propeller geometry and inflow nominal wake (Choi & Kinnas 1998a, Kinnas et al 1999). In that method, the
propeller effect, which appears through the body force terms in the Euler equations, is averaged with time. Nevertheless,
the propeller loading is allowed to vary in space with strength which depends on the loading at each blade angle.
In the present work, a fully three-dimensional unsteady Euler solver, based on a finite volume approach and the
pressure correction method, is developed and applied to the prediction of the unsteady effective wake for propellers
subject to non-axisymmetric inflows. The proposed unsteady analysis, which is the unsteady extension of our previous
work (Choi & Kinnas 1998a), captures the truly unsteady behavior of the interaction between the vorticity in the inflow
and the propeller action by adopting the body force distribution which varies in both space and time.
2 PRESENT METHOD
The unsteady flow field around a propeller can be decomposed into two parts; one rotational and one irrotational
field. In the present method, the propeller induced velocity field, is the irrotational velocity field. The propeller
induced velocity, can be expressed in terms of the perturbation potential, by the following relation.
(1)
The propeller vortex lattice method (VLM), named MPUF-3A, is used to solve for the perturbation potential on the
propeller as shown on the left side of Figure 1 (Kinnas et al 1998, Kinnas et al 1999).
the authoritative version for attribution.
OCR for page 617
lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
AN UNSTEADY THREE-DIMENSIONAL EULER SOLVER COUPLED WITH A CAVITATING PROPELLER ANALYSIS METHOD 617
The perturbation potential as defined by equation (1) with respect to a coordinate system which rotates with the blade at a
constant angular velocity, should satisfy the Laplace equation. Note that the propeller rotation, has only an axial
component which is negative for right-handed propellers. The boundary condition on the blade mean camber surface is
Figure 1: Interaction between the potential flow vortex lattice method (VLM) and the Euler equation based finite volume
method (FVM).
(2)
where, is the unit vector normal to the mean camber surface, and is the effective wake as defined in the next
paragraph. In addition, a Kutta condition must be applied at the trailing edge of the blades with the assumed trailing wake
geometry.
Then, the rotational part can be expressed by subtracting this irrotational field, from the total velocity field,
The rotational part determined in this way is defined as the effective wake velocity, given as follows.
(3)
The effective wake, is used as the inflow, when solving for the propeller perturbation potential, The total
velocity, is evaluated from the Euler solver, as described in the next paragraphs. Once the potential is found, two
things can be calculated; (a) the propeller loading as a function of space and time, and (b) the propeller induced velocity
or equivalently The propeller loading is needed to determine the body force distribution in the Euler solver, and
the propeller induced velocity is needed when the effective wake is calculated. If unsteady sheet cavitation exists on the
blade, then the extent and the thickness of the cavities are also determined as part of the solution.
On the other hand, the global flow field is rotational and is solved by the finite volume method Euler solver, as
shown on the right side of Figure 1. The propeller loading obtained from the VLM is converted into the body force and
this body force is distributed over the cells corresponding to the blade location. When the Euler equations,
(4)
are solved, the total velocity is obtained. The effective wake can then be computed by equation (3) using the
propeller induced velocity, already known from the propeller VLM. The overall process is iterative, between the
potential flow solution of the propeller and the Euler equation solution of the global flow field, as shown in Figure 1.
The main assumption behind the current approach is that viscosity plays a less significant role in the interaction
between the rotational inflow and the propeller action, even though it plays a major role in developing the nominal wake.
This is a valid assumption in the vicinity of the propeller where the propeller-inflow interaction is dominated mainly by
inviscid vorticity dynamics. In order to predict the nominal wake, a high Reynolds number viscous flow around the ship
hull must be solved, and this sometimes has been found to be a very difficult task (Larsson 1997). Instead, one can use the
measured nominal wake so that the “real fluid effects” (including the effects of viscosity) are included in the nominal
wake data. The interaction between the vorticity in the nominal wake and the action of the propeller can then be captured
by solving the Euler equations, and the effective wake can be calculated from the difference between the total flow field
and the flow field induced by the propeller.
In the current algorithm, all components of the inflow must be specified at the upstream boundary of the domain in
which the Euler equations are solved. This boundary must be placed at a sufficient distance upstream of the propeller
plane so that the specified inflow is not affected by the presence of the propeller. Ideally, the nominal wake inflow (in the
absence of the propeller) must be measured at the upstream boundary of the calculation. However, often the nominal
wake is measured at the propeller plane, and thus the inflow at the upstream boundary needs to be deduced from the
the authoritative version for attribution.
measured. This situation was also encountered in the effective wake analysis of (Huang & Groves 1980). In some cases,
OCR for page 618
lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
AN UNSTEADY THREE-DIMENSIONAL EULER SOLVER COUPLED WITH A CAVITATING PROPELLER ANALYSIS METHOD 618
a reverse flow calculation is attempted in order to find the flow at the upstream boundary which will produce the
measured nominal wake at the propeller plane. It is concluded that placing the nominal wake at the upstream boundary is
the easiest approach, especially when the vertical component of the inflow is small compared to the axial component. One
more variation in the current method is that the effective wake is evaluated not at the propeller plane but at a plane just
upstream of the propeller, due to the difficulty in computing the induced velocity near the propeller panels. Luckily, the
location where the effective wake is evaluated has been found to have insignificant effect on the prediction (Choi &
Kinnas 1998a, Choi 2000).
Figure 2: Ship-fixed Cartesian coordinate system.
3 FORMULATION
A ship-fixed Cartesian coordinate system, as shown in Figure 2, is used in the present formulation. The origin of the
coordinate system is at the center of the propeller disk, with the x-axis pointing downstream along the propeller shaft axis.
In the coordinate system, the positive y points vertically upward, and the positive z-axis port side (left) of the ship. Also
the angle θ is measured around the x-axis starting from the y-axis, as shown in Figure 2.
The variables are made non-dimensional by a reference length which is chosen as the propeller radius R, and a
reference velocity equal to the velocity at infinity U∞ (ship speed).
(5)
(6)
(7)
Figure 3: Cell (i, j, k) on which the finite volume method (FVM) is applied.
(8)
where, the hat(ˆ) means a dimensional variable, and is the body force per unit mass.
3.1 Steady Euler Solver
In the steady Euler solver, named WAKEFF-3D, the method of artificial compressibility (Chorin 1967) is applied. In
the artificial compressibility method, pseudo-unsteady terms are introduced on the left hand side of the steady
incompressible Euler equation to mimic the solution procedure of the unsteady compressible Euler equations.
After introducing the pseudo-unsteady terms, the Euler equations can be written in the following three-dimensional
Cartesian form using nondimensional variables.
(9)
where, t* is a pseudo-time which can be regarded just as an iteration parameter. The terms U, F, G, H, and Q are
the authoritative version for attribution.
defined as follows.
(10)
where, β is the artificial compressibility factor which is a controllable constant. If steady state is reached, the first
term in equation (9) vanishes and the incompressible steady flow solution is obtained. For the rest of this section, the
superscript, *, on the time
OCR for page 619
lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
AN UNSTEADY THREE-DIMENSIONAL EULER SOLVER COUPLED WITH A CAVITATING PROPELLER ANALYSIS METHOD 619
variable, t, will be omitted with the understanding that it is the pseudo-time.
Figure 4: Arrangement of the nodes and cells in the present method.
According to the finite volume method (FVM), the Euler equations (9) are integrated over a cell of volume V. Then,
with the use of the divergence theorem, the integral form of the Euler equations is applied to each cells as shown in
Figure 3. The semi-discrete equation in which only the space is discretized can be written as follows for each cell.
(11)
where, Vijk is the volume of the cell (i, j, k). The areas (Ax, Ay, Az)f are the projections of the area of the face, f, along
the (x, y, z) directions, respectively. Note that the vertex based scheme, in which all variables in U are stored at the nodes,
is used in the present method.
For the time discretization, Ni's Lax-Wendroff method (Ni 1982) is applied. That is, the variable U at node, “5”
(shown in Figure 4), at the next time (pseudo-time) step, n+1, is approximated by the following second order formula.
(12)
where, ∆t is the time step size, and the superscript n represents the value of the current time step. The first order
term, ∂U/∂t, and the second order term, ∂2U/∂t2, are approximated by the average of their values at the neighboring cells,
where the values at each cell are expressed using the discretized governing equation (11) itself. Then, the resultant
equation becomes the distribution formula, in which the contributions from each neighboring cell are grouped together.
The contributions from the surrounding cells to the node “5”, are given in (Choi 2000) with more details of the
formulation.
The artificial dissipation (or viscosity) is applied to improve the stability of the numerical method (Anderson 1995).
The second and fourth order dissipations, respectively µ 2 and µ 4, are scaled by ∆t and added to the distribution formula.
Figure 5: Finite volume method applied to the cylindrical grid in the Cartesian coordinate system.
(13)
with
(14)
(15)
where, δii, etc. are the finite central difference operators. For example, the fourth order i-directional difference
operator, δiiii, is defined as follows.
(16)
The artificial dissipation coefficients, σ2 and σ4, are user specified constant parameters that control the amount of the
dissipation. For most of the cases presented in this paper, the second order dissipation coefficient, σ2, is zero. The time
step size should be determined so that the CFL condition (Courant et al 1967) is satisfied. The iteration in time (pseudo-
time) is continued until the maximum change of variables is less than a certain tolerance.
Since a cylindrical grid system fits the propeller problem naturally, a cylindrical grid system is used in discretizing
the authoritative version for attribution.
the fluid domain, as shown in Figure 5, while a Cartesian coordinate system is used to describe the space variables and the
velocity vectors, as well as the Euler equations. The k-index of the grid is increasing in the θ-direction and the surface
corresponding to the last k-index and that of the first k-index are the same surface. On this surface, a periodic boundary
condition is applied.
3.2 Unsteady Euler Solver
In the unsteady Euler solver, named WAKEFF-3U, the finite volume method is applied only for the dis
OCR for page 620
lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
AN UNSTEADY THREE-DIMENSIONAL EULER SOLVER COUPLED WITH A CAVITATING PROPELLER ANALYSIS METHOD 620
cretization of the three unsteady momentum equations, while the continuity equation is handled by the pressure correction
method (Patankar 1980, Rhie & Chow 1983).
3.2.1 Finite Volume Discretization for Momentum Equations
In the unsteady formulation, the following dimensionless incompressible Euler (momentum) equations in differential
form are used.
(17)
where, the column matrices U, F, G, H and Q, each of which has only three elements, are defined as follows.
(18)
Note that the pressure gradient terms are grouped in Q with the body force terms because the pressure correction
method is adopted in the unsteady formulation. From this point, the formulation and the numerical implementation of
equation (17) are very similar to those of the steady case, described in the previous section, and are not repeated here.
3.2.2 Pressure Correction Method for Continuity Equation
In the present method, the idea of the SIMPLE method (Patankar 1980) is followed. As opposed to the staggered grid
system of the original SIMPLE method, the collocation scheme, in which all variables are defined and stored at nodes, is
applied. To achieve this, the pressure weighted interpolation (Rhie & Chow 1983) is applied when the pressure gradient is
discretized.
Using equation (13) and the definition of U by equation (18), the velocity at the next time step at node “5” can be
expressed as follows.
(19)
Note that the artificial dissipation is based on the velocities at time step n, and that the terms include the
residual, Rc, and the source, Qc, of the surrounding cells.
The intermediate velocity field denoted by the superscript, *, is defined as the one based on an intermediate (or
guessed at first) pressure field, p*. Since the method is iterative, to find the correct (within a certain tolerance) velocity
and pressure fields for time step n+1, the intermediate fields are updated repeatedly. Thus, the intermediate velocity and
pressure are not marked with the time step index, understanding that they are between the time steps n and n+1.
(20)
where, is evaluated by using the intermediate pressure field, p*. Also, the artificial dissipations, and
are evaluated based on the intermediate velocity field, U*.
The pressure and velocity corrections, which are represented by prime (′), are defined as follows.
(21)
(22)
where, and are the corrections to the pressure and the velocity fields, respectively, so that satisfies the
continuity equation. According to the derivations given in (Choi 2000), the velocity correction, is expressed by the
pressure correction and the pressure correction can be found by solving a discretized equation for the pressure
correction. The pressure correction equation can be solved iteratively using the following explicit expression for with
the under-relaxation, αpp.
(23)
The optimum value of the under-relaxation, αpp, is found to be 0.5 in the present method. In summary, the following
operations are performed within each time step of the present unsteady method.
Sequence of Operations
the authoritative version for attribution.
1. Guess the pressure field p*.
2. Solve the momentum equation by equation (20), to obtain intermediate velocities u*, υ* and w*.
3. Solve for the pressure correction, p′, by iteratively applying the equation (23). The optimum number of this
inner iteration is found later to be IppIter=2.
OCR for page 621
lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
AN UNSTEADY THREE-DIMENSIONAL EULER SOLVER COUPLED WITH A CAVITATING PROPELLER ANALYSIS METHOD 621
4. Calculate the (correct) pressure pn+1 using equation (21), and the (correct) velocities un+1, υn+1 and wn+1 using
equations (22).
5. Treat the corrected pressure pn+1 as a new intermediate pressure, p*, return to step 2, and repeat until either
the residual of the continuity equation, QM, is smaller than a given tolerance or a prescribed number of
iterations, IpcMax, is reached.
Later, the optimum number of pressure correction per each time step is determined to be IpcMax=10.
3.3 Boundary Conditions
The six boundary conditions applied in the present methods are summarized next.
• Upstream boundary:
(24)
• Downstream boundary:
(25)
• Axisymmetric shaft or hub boundary:
The free slip condition is applied by using the tangential velocity component, uθ.
(26)
where, the unit normal vector, (nx, nr), pointing into the fluid region is defined to be positive.
• Center line boundary (j=1 line in the grid):
(27)
where, Nk is the number of nodes in the circumferential direction.
• Outer boundary (far field):
(28)
• Periodic boundary:
(29)
where, the k-index is shown in Figure 5.
Figure 6: A cell of the finite volume method (FVM) mesh that contains the body surface.
4 BODY FORCE
To find an appropriate distribution of body force which represents the propeller (or hydrofoil in two-dimensions) in
the Euler solver, two-dimensional lifting body problems are discussed first.
4.1 Two-Dimensional Body Force
For simplicity a steady incompressible flow around a two-dimensional hydrofoil is considered for the moment.
Figure 6 shows the part of the body surface intersected by a cell of the FVM mesh. Pressure p is acting on the part of the
body surface of length λ, and the cell has a size of ∆x×∆y. From the dimensionless Euler equations for two-dimensional,
incompressible, steady flows, the following expression for the body force can be derived for the cell shown in Figure 6.
(30)
(31)
where, the subscripts N, S, W, and E stand for the four edges of the cell at north, south, west, and east, respectively.
Note that the first two terms of the discretized equations (30) and (31) are the forces due to the momentum convection.
The last terms of the discretized equations above can be interpreted that the body force at a cell should be equivalent
to the integrated pressure force acting in the normal direction over the part of the body surface contained in the same cell.
the authoritative version for attribution.
When the cell shown in Figure 6 is assumed to be sufficiently small, the body force can be expressed as follows by
equating the two forces; the total force due to the pressure, and the integration of the body force over the cell,
(32)
If the body surface is vertical, the body force simply becomes fx=p/∆x, and if the body surface is horizontal, the body
force becomes fy=p/∆y. These
OCR for page 622
lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
AN UNSTEADY THREE-DIMENSIONAL EULER SOLVER COUPLED WITH A CAVITATING PROPELLER ANALYSIS METHOD 622
expressions for the (x, y)-directional body forces are consistent with the last terms of the discretized equations (30) and
(31) under the condition that the pressure inside the body is assumed to be zero.
Figure 7: Stream lines and pressure contours of the potential flow solution; NACA0015 thickness, NACA a=0.8 camberline,
5% maximum camber, 5° angle of attack.
Figure 8: Body force cells and the body force vectors that represent the foil according to the surface distribution model.
(Arrows represent the force acting on the fluid at each of 292 body force cells.)
4.1.1 Surface Distribution Model
In the surface distribution model, the body force is applied to the cells that contain a part of the boundary element
method (BEM) body surface element. The first and second terms of the equations (30) and (31) are calculated by using
the local FVM velocity (u, v), so that they can exactly cancel the convective terms on the left hand side of the momentum
equations, when the body forces are substituted into the momentum equations. On the other hand, the last terms are
computed from the surface pressure known from the BEM solution by applying equation (32) and by setting zero pressure
at the nodes inside the body.
Figure 9: Stream lines and pressure contours of the Euler solution (surface distribution model); NACA0015 thickness, NACA
a=0.8 camberline, 5% maximum camber, 5° angle of attack.
(33)
When this body force model is applied, the momentum equations provide a direct balance between the surface
pressure force and the cell body force. In other words, this method is equivalent to specifying the pressure of the cell as
the authoritative version for attribution.
the value obtained from the boundary element method (BEM) solution at the same location.
The surface distribution model is applied to a two-dimensional foil made of a NACA0015 thickness and a NACA
a=0.8 camberline with 5% maximum camber. The geometry of the foil is represented by 100 boundary elements in the
BEM solution. In Figure 7, the potential flow solution around the foil is shown for uniform inflow with an angle of attack
of 5°. The potential flow solution is obtained from the two-dimensional perturbation potential method described and fully
verified in (Choi & Kinnas 1998b). The flow field very near the surface of the foil is not accurate because the evaluation
of the induced velocity at locations which are closer to the BEM element
OCR for page 623
lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
AN UNSTEADY THREE-DIMENSIONAL EULER SOLVER COUPLED WITH A CAVITATING PROPELLER ANALYSIS METHOD 623
than the size of the element, is inaccurate due to the constant strength singularity approximation. However, note that the
pressures on the surface control points which are used in the calculation of the FVM body force are accurate because they
are obtained directly, by differentiating the solution from the BEM along the surface of the body.
In Figure 8, the 292 body force cells and the body force vectors at each cell, that represent the foil, are shown. It is
observed in the figure that the largest part of the suction side is pulling the fluid while the largest part of the pressure side
is pushing the fluid.
In Figure 9, the solution of the Euler equations is shown, as obtained with the body force distribution shown in
Figure 8. The formulation of the two-dimensional Euler solver is similar to that of the three-dimensional problem
explained earlier, when the z coordinate is omitted. The definition of the artificial dissipation is slightly different only for
this specific case. That is, the last term in equation (13) should be replaced with ∆t(µ 2−µ 4)/Sij, where Sij is the average of
the areas of the cells neighboring the node (i, j). For this computation, the square domain extending −8≤x≤8 and −8≤y≤8
is used with the boundary conditions specifying (u, v, p) at the upstream, top, and bottom boundaries and ∂( u, v, p)/∂x=0
at the downstream boundary. The pressure field and the streamlines obtained from the Euler solver are very close to those
of the potential flow solution shown in Figure 7. For the most part of the field, excluding regions very near the foil
surface and the wake surface, the difference between the two velocity fields is observed to be less than 1% of the inflow
velocity. It is an important fact that the flow field, including the effect of the foil thickness, is well simulated with the
surface distribution model. Also, note that the FVM computation required about 5 hours of CPU time on a DEC Alpha
Station 600 5/2661, while the BEM computation took only 0.1 second on the same machine.
4.1.2 Camberline Pressure Model
It is natural that very thin foils should be represented by the body force cells along the mean camber surface of the
foil. This can be achieved easily by just relocating the pressure gradient part of the body force determined by equation
(32) for the cells on the BEM surface to the cells on the mean camber surface, resulting into two cell thickness
distribution of body forces. A simpler version of this method is to have only one cell thickness distribution along the
camberline by adding the two force vectors from the upper and the lower surface of the foil. Numerical tests showed that
both body force cell arrangements produced practically the same velocity fields with only small local differences (Choi
2000). The one cell thickness distribution has the advantage that there is no need to set zero pressure inside the foil
because the addition of the upper and the lower pressures automatically implies the assumption that the two pressures are
based on a certain common reference value, say zero, and there is no need to set this explicitly. It can also be shown that
the finite volume discretization of a cell that contains both surfaces of the foil arrives at the same conclusion.
The above mentioned body force model can be made even simpler by neglecting the convective terms which used to
be calculated from the local FVM solution, and keeping only the pressure gradient term, as shown next.
(34)
where, the superscript + and − mean the values from the upper surface and the lower surface of the BEM boundary,
respectively. This body force model which is simpler to apply is named as the camberline pressure model. It has been
confirmed by numerical tests that the original camberline distribution model (with both FVM and BEM terms) and the
simpler camberline pressure model (with only BEM terms) do not make any significant difference in the predicted
velocity at the locations corresponding to the effective wake evaluation plane in propeller applications.
To test the camberline pressure model, a relatively thin hydrofoil geometry which would be similar to a section of a
propeller at outer radii is selected. A two-dimensional foil section constructed from NACA0005 thickness and NACA
a=0.8 camberline with 2% maximum camber is used in this test.
The potential flow BEM solution around the foil, which is represented by 100 elements, is shown in Figure 10 when
uniform inflow with 5° angle of attack is incident. The potential flow solution is obtained again from the two-dimensional
perturbation potential method described in (Choi & Kinnas 1998b). Again, the flow field very near the surface of the foil
is not accurate because of the inaccuracy of the induced velocity calculation routine.
The body force distribution is shown in Figure 11, where only 10 body force cells are used to represent the foil. Note
that the starting point of each force vector is the lower left corner of the cell which the body force belongs to. The body
force has smaller magnitude than that from Figure 8 because what
the authoritative version for attribution.
1DEC Alpha Station 600 5/266 is approximately equivalent to a 500 MHz Pentium PC.
OCR for page 624
lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
AN UNSTEADY THREE-DIMENSIONAL EULER SOLVER COUPLED WITH A CAVITATING PROPELLER ANALYSIS METHOD 624
is plotted is the body force per unit mass and the cell area where the force is to be distributed becomes larger now.
Figure 10: Stream lines and pressure contours of the potential flow solution; NACA0005 thickness, NACA a=0.8 camberline,
2% maximum camber, 5° angle of attack.
Figure 11: Body force cells and the body force vectors that represent the foil according to the camberline pressure model,
equation (34); grid with 50x50 cells (10 body force cells).
In Figure 12, the solution of the Euler equations is shown, which is obtained with the body force distribution shown
in Figure 11. The formulation of the two-dimensional Euler solver is similar to that of the three-dimensional problem
including the way the artificial dissipation is defined. The domain size and the boundary conditions are the same to those
applied in the previous case. A relatively coarse grid with 50x50 cells is used because the number of body force cells in
this grid is similar to the number of cells along the chord of a blade in the propeller flow problems. The jagged pressure
contours near the foil could be smoothed out by using increased artificial dissipation. However, the present result with the
moderate artificial dissipation is preferred because the present artificial dissipation seems to be adequate in the region of
interest, i.e. not too close to the camberline. The pressure field and the streamlines obtained from the Euler solver are
similar to those of the potential flow solution shown in Figure 10, even though the foil thickness effect is completely lost
in the Euler solution.
Figure 12: Stream lines and pressure contours of the Euler solution (camberline pressure model); NACA0005 thickness,
NACA a=0.8 camberline, 2% maximum camber, 5° angle of attack.
The difference of the flow fields computed by the potential solver and the Euler solver is shown in Figure 13. The
difference of the two solutions is computed by subtracting the potential flow field (Figure 10) from the Euler solution
flow field (Figure 12). For the most part of the field excluding very near the foil surface and the wake surface, the
difference between the two solutions is less than 1% of the inflow velocity.
The calculation of the effective wake requires the evaluation of velocities on a plane at about 0.3 propeller radius
upstream of the propeller plane. If a section of a propeller blade at around 70% of the propeller radius is considered, such
plane would correspond to a straight line as the dashed line shown in Figure 13, which is sloped relative to the nose-tail
line of the foil section. Along this line, the velocity magnitudes of the Euler solution and the potential solution are
the authoritative version for attribution.
compared in Figure 14. Along this line, the maximum change of the velocity is about 7% of the inflow velocity, and the
maximum error of the Euler solution is about 1%. The computing time with this 50x50 cell mesh is less than 3 minutes on
a DEC Alpha Station 600 5/266.
OCR for page 625
lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
AN UNSTEADY THREE-DIMENSIONAL EULER SOLVER COUPLED WITH A CAVITATING PROPELLER ANALYSIS METHOD 625
Figure 13: Difference of the velocity magnitude between the Euler solution (camberline pressure model) and the potential
solution, i.e. (Euler solution—potential solution); The dashed line from (−2, 0) to (1, 1.5) is the cut along which the two
solutions and the difference are compared in Figure 14.
Compared to the run on a 200x200 mesh with surface distribution model presented in the previous section, the
current run on a 50x50 mesh with the camberline pressure model requires only 1/100 of the computing time, while the
predicted velocity field at locations of interest has the same 1% maximum error. A three-dimensional version of the
camberline pressure model with relatively coarse mesh is applied to compute the body force in the unsteady propeller
problem.
4.2 Body Force Representation of an Actuator Disk
In actuator disk theory (Breslin & Andersen 1994), the (dimensional) thrust, T, is assumed to be uniform over the
disk area, Ao. The thrust loading coefficient, CT, which is based on the ship speed, U∞, is defined as follows.
(35)
If the body force is distributed over the dimensionless distance ∆x in the axial direction, the dimensionless body
force per unit mass, fx, is expressed as follows.
(36)
When the thrust loading coefficient, CT, is specified, the body force, fx, can be readily computed by equation (36).
The axial distance ∆x is chosen to be two cell thickness, after it has been confirmed that there is only local difference in
the solution if this distance is varied from one cell to four cell thickness for the usual mesh of which the cell size near the
disk is 0.04.
Figure 14: Comparison of the velocity magnitudes predicted by the Euler solver and the panel method along the dashed line
shown in Figure 13; 50x50 cells, body force from the camberline pressure model, using equation (34).
4.3 Steady Body Force Model
In the steady three-dimensional Euler solver, the body force which represents the propeller can be a function of
space, i.e. The variation in the axial direction is neglected and the body force is distributed over the disk area of
the propeller within two cell thickness axially. The body force distributed in this way can simulate the time averaged
spatially varying loading, such as in the case of a propeller subject to inclined uniform inflow.
The body force distribution over the disk area is calculated in the lifting line sense. First, the instantaneous blade
circulation distributions Γ(r) along the lifting lines at each blade angle are computed by the cavitating potential flow
solver (Kinnas et al 1998). Then, the circulation distribution Γ(r, θ) over the disk area is obtained by superimposing all of
the instantaneous blade circulation distributions Γ(r), considering the skew. Finally, it is assumed that the radial and
circumferential distribution of the thrust and the torque follow the radial and circumferential distribution of the circulation Γ
the authoritative version for attribution.
(r, θ). In this work, only the axial and tangential body forces, fx and fθ are applied. The details of the steady body force
model can be found in (Choi & Kinnas 2001).
OCR for page 626
lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
AN UNSTEADY THREE-DIMENSIONAL EULER SOLVER COUPLED WITH A CAVITATING PROPELLER ANALYSIS METHOD 626
Figure 15: Unsteady body force model for a propeller blade; the cylindrical surface at a radius is expanded; a view from top.
4.4 Unsteady Body Force Model
In the unsteady solution, the body force must be a function of time as well as space. Since the blade-to-blade flow
field in the vicinity of the propeller is resolved in the unsteady solution, the blade should be represented with more
accuracy than that used in the steady solution.
Three unsteady body force models were considered in this work; (a) the body force distributed over the flat propeller-
plane-projected blade area, proportional to the blade bound vortex distribution Γ(r, θ, t), (b) the body force model utilizing
the bound vortex distribution on blade which includes the axial variation, Γ(x, r, θ, t), and (c) the unsteady three-
dimensional extension of the camberline pressure model explained in section 4.1.2. The model (c) was finally selected
because numerical studies showed that the models (a) and (b) gave unsatisfactory total velocity (about ±5% error) at the
plane where the effective wake is evaluated. The selected model has the advantage that it is based on the blade pressure
distribution which is directly related to the blade force and, thus, more rigorous (also more accurate) than the other two
models based on the blade circulation distribution. Note that a time averaged version of this method is conceptually
similar to the steady body force model of (Toda et al 1990).
The unsteady body force for propeller blades is obtained by an extension of the two-dimensional camberline pressure
model. The pressure, the unit normal vector, and the area of the BEM element are known for each side of the blade from
the propeller potential flow solution. The product of the pressure, the unit normal vector, and the area is calculated first
for each side, and then the two force vectors from each side are added to obtain the force vector corresponding to the
element on the camber surface. This is done for the 60 blade angles spaced by 6°, so that the pressure force and the
camber surface geometry can be found by interpolation from these 60 sets of data for the blade angle corresponding to the
non-dimensional time clicking in the Euler solver. After the time-interpolation of the pressure force has been completed,
the camber surface geometry and the pressure force are interpolated in the radial direction first. For a certain radius, the
camber surface appears as a camberline overlapped by the Euler solver mesh as shown in Figure 15. Then, the body force
for the cells for which the camberline goes through can be computed by collecting the pressure force over the segment of
the camberline contained by each cell.
5 NUMERICAL RESULTS
Several convergence studies were performed to find the effect of the grid resolution and the effect of the artificial
dissipation coefficient, as reported in (Choi 2000, Choi & Kinnas 2001). A grid with 80 cells in the axial direction and 50
cells in the radial direction, covering from four propeller radii upstream to four propeller radii downstream and expanding
five propeller radii radially, was found to be appropriate. Also, zero second order artificial dissipation and fourth order
artificial dissipation coefficients between 0.5 and 1.5 were found to be appropriate for the examined cases. Concerning
the circumferential number of cells, it was found that at least 12 cells were required to resolve each harmonic of the
effective wake; i.e. 60 cells were required to resolve 5 harmonics.
5.1 Actuator Disk
An actuator disk with thrust loading coefficient CT=3.0 is analyzed using the unsteady Euler solver. In uniform flow,
the uniform axial body force is applied suddenly on the cells at the disk at the initial dimensionless time t=0 and the time
evolution of the solution is examined. The final solution reached at the end of the transient flow field should agree with
the steady actuator disk result.
To solver this unsteady actuator disk problem, the grid with 80x50x30 cells is used with the following parameters;
5000 time steps with ∆t=0.002, σ2=0, and σ4=1. The run takes 26.4 CPU hours on a DEC Alpha Station 600 5/266, and
the convergence is shown in Figure 16. The curves labeled as ∆u etc. stand for the maximum (over all cells) value of the
change between time steps of the axial velocity u etc. at each node. The curve labeled as QM in this figure represents the
maximum error in satisfying the continuity equation, i.e. the residual of the continuity equation. The dimensionless time,
the authoritative version for attribution.
t, is scaled
OCR for page 628
lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
AN UNSTEADY THREE-DIMENSIONAL EULER SOLVER COUPLED WITH A CAVITATING PROPELLER ANALYSIS METHOD 628
be observed that the axial velocity curve reaches its final shape at around t=4. It can be observed from the figure that the
converged pressure and velocity along the center line agree well with the analytic values.
Figure 19: Time evolution of the z-vorticity field; t=1.0 (top), t=2.0 (bottom).
The time evolution of the solution field on the center plane is shown in Figure 18. It is observed that the entire flow
field has reached almost a steady state near t=4. In Figure 19, the z-vorticity (ωz) fields are shown on the center plane
corresponding to the dimensionless times t=1.0 and t=2.0. It is clearly shown in the figure that the start-up vortex ring
which originates from the rim of the actuator disk at t=0.0, evolves forming a contracting vortex tube. Since the axial
velocity inside the slip stream is faster than that outside, the z-vorticity, ωz, has positive value in the upper half (y>0) of
the domain, and negative in the lower half of the domain, as can be seen in the figure. As the start-up vortex convects
downstream, it diffuses due to the coarser cells downstream and due to the artificial dissipation.
There are two kind of iterations in the present unsteady method; one is the pressure correction iteration, of which the
maximum number is IpcMax, and the other is the inner iteration to solve the pressure correction equation of which the
maximum number is IppIter. To find the optimum number of iterations, several combinations of IpcMax and IppIter are
tried on the 50x30x30 grid, and the CPU times for each run are compared in Table 1. It is observed from Table 1 that too
many IppIter can cause the run to diverge. The solution and the convergence history of the converged runs are found to
have little dependence on IppIter. On the other hand, IpcMax is found to have direct effects on the level of QM, though the
gain in accuracy becomes very small for IpcMax≥10. Compromising between the run time and the convergence, the
optimum values have been chosen as IppIter=2 and IpcMax=10.
Table 1: CPU hours on a DEC Alpha Station 600 5/266 for combinations of IpcMax and IppIter.
IpcMax IppIter
2 4 6 8 12 18
2 1.0 1.0 1.1 × ×
4 1.9
6 2.7 2.9 3.2 × × ×
10 4.5 × × ×
12 5.9
18 8.0 8.9 9.5 × ×
20 9.0
(× means the run diverged.)
5.2 Uniform Inflow
Due to the fact that the body force cells are placed along the blade camber surface rather than on the “exact” blade
the authoritative version for attribution.
surface, the blade thickness effect on the flow field is completely lost in the unsteady Euler solver, as observed in the two-
dimensional tests. This means that the flow field obtained by the unsteady solver is not consistent with the flow field
predicted by the propeller potential flow solver simply because the propeller potential flow solver does include the blade
thickness effect. Placing the body force cells on the “exact” blade surface, using the present discretization in the Euler
solver, is practically impossible due to the fact that the resulting size of a cell is comparable to the blade thickness.
Therefore, it is necessary to quantify this inconsistency because the difference between the two flow fields is directly
related to the accuracy of the predicted effective wake.
The same inconsistency should affect the steady Euler solver since it uses the disk-like distribution of the body force
based on the blade circulation distri
OCR for page 629
OCR for page 631
OCR for page 632
OCR for page 633
OCR for page 634
OCR for page 635
OCR for page 636
OCR for page 637
OCR for page 638
lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
AN UNSTEADY THREE-DIMENSIONAL EULER SOLVER COUPLED WITH A CAVITATING PROPELLER ANALYSIS METHOD 629
bution. However, even though the blade thickness effect cannot be properly treated by the steady body force model, the
validation test for the uniform inflow, presented in (Choi & Kinnas 2001), shows that the difference between the two flow
fields, predicted by the steady Euler solver and the potential flow solver, is as small as 1% of the inflow velocity. This
small error can be explained from the fact that the time averaged effect of the blade thickness is small as compared to the
local instantaneous effect which can be larger. In other words, the blade thickness effect becomes more of an issue only in
the unsteady Euler solver because the instantaneous blade-to-blade flow field, which is affected by the thickness of blades
nearby, is determined in the unsteady Euler solver.
Figure 20: Axial velocities predicted by the unsteady Euler solver and by the potential flow solver for the one-bladed DTMB
4148 propeller; xeff =−0.323.
In this section, only the effect of the blade thickness on the blade itself is examined by considering an one-bladed
propeller, so that the influence from the thickness of the other blades of the propeller is avoided. The one-bladed propeller
is made by removing two of the three blades from the DTMB 4148 propeller, while keeping the geometry of the
remaining blade (including the thickness) the same. The geometry of DTMB 4148 propeller can be found in (Choi 2000).
The inflow condition is unit uniform inflow so that the velocity field predicted by the potential theory can be compared
directly with that predicted by the Euler solver. The propeller is working at the advance ratio, Js≡nD/U∞=0.8, under the
fully wetted condition, where n is the number of propeller revolution per second and D is the propeller diameter. The grid
has 80x50x60 cells in the domain −4
lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
AN UNSTEADY THREE-DIMENSIONAL EULER SOLVER COUPLED WITH A CAVITATING PROPELLER ANALYSIS METHOD 630
locities along the same circumference are compared in Figure 22. The agreement between the two velocities is very good.
As a result, the effective wake has a smaller error ranging from −0.7% to 0.2%. When the MPUF-3A results in the two
Figures 20 and 22 are compared to each other, the blade thickness effect is more prominent in the deceleration region
(near θ=−25°) than in the acceleration region. It is very interesting that the maximum error of the original thickness result
shown in Figure 20 occurs also near the deceleration region. In the deceleration region, the velocities predicted by the
Euler solver are found to be affected very little by the blade thickness, while the potential solutions show appreciable
effect of the blade thickness. Therefore, it can be concluded that the disagreement near the deceleration region shown in
Figure 20 is mostly due to the defect of the blade thickness effect in the Euler solver.
Figure 22: Axial velocities predicted by the unsteady Euler solver and by the potential flow solver for zero-thickness one-
bladed DTMB 4148 propeller; xeff=−0.323.
In Figure 23, the velocity fields on the plane at 0.323 radius upstream of the propeller are shown for the zero-
thickness propeller case. It is observed that the total and the propeller induced velocities have very similar contour shapes,
and that the effective wake has maximum error near the hub ranging ±1%. The error near the hub is contributed by the
difference of boundary condition at the hub. The zero normal velocity boundary condition is applied on the hub in the
Euler solver, while the hub is not included in the potential flow solver. Even though the potential flow solver has the
ability to consider the hub via a simplified image model, the hub is turned off in the present potential solver computation
because the simplified image model has not been found to produce the correct flow field upstream of the propeller plane.
Figure 23: Comparison among the total, the propeller induced, and the effective wake as predicted by the unsteady Euler
solver at 0.323 radius upstream of the zero-thickness one-bladed DTMB 4148 propeller.
To confirm the blade thickness is responsible for the inconsistency shown in Figure 20 again, the following test is
performed: When the propeller induced velocity is computed, the velocity induced by the dipoles (equivalent to vortex
loops) and sources on the blades and the associated trailing wakes are summed. The dipoles are responsible for the lift of
the blade, while the sources represent the blade thickness and the cavity thickness, if any. If the velocity influence due to
the sources is omitted in the summation, the blade thickness effect on the induced velocity is excluded, and thus the
potential flow model becomes consistent to the camberline pressure model in the Euler solver. The axial velocities
predicted in this way are shown in Figure 24. Note that the agreement between the Euler solution and the potential flow
solution is as good as in the case of zero-thickness propeller, shown in Figure 22.
In conclusion, the effect of the blade thickness is found to be the major source of error in the calculation of the
effective wake. Once the blade thickness is eliminated either by using a zero-thickness propeller or by omitting the
velocity influence of the thickness sources, the error in the effective wake is less than 1% of the inflow velocity.
5.3 Propeller DTMB 4842 in Uniform Inflow
For the unsteady analysis of the five bladed DTMB 4842 propeller, working at Js=0.9 in uniform inflow, a grid with
80x50x120 cells which covers the domain −4≤x≤4, 0.323≤r≤5 is used. The propeller
the authoritative version for attribution.
lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
AN UNSTEADY THREE-DIMENSIONAL EULER SOLVER COUPLED WITH A CAVITATING PROPELLER ANALYSIS METHOD 631
geometry is given in (Choi 2000). The number of cells in the circumferential direction is chosen to be 120 because this
propeller has five blades and more variation is expected in the circumferential direction. The blade thickness is reset to
zero. One revolution of the propeller is divided into 900 time steps using a time step size ∆t=0.002, the unsteady
computation parameters are IpcMax=10, IppIter=2, αpp=0.5, and the artificial dissipation coefficients are σ2=0 and σ4=5.
The numerical simulation for the one revolution of the propeller took 39.5 CPU hours on a DEC Personal Workstation
500au2.
Figure 24: Axial velocities predicted by the unsteady Euler solver and by the potential flow solver for the one-bladed DTMB
4148 propeller at xeff =−0.323. When the propeller induced velocity is computed, the influence of the blade thickness is
excluded.
In Figure 25, the axial body force distribution that represents the propeller in the unsteady Euler solver is shown at
the instance the key blade is at the top position (0º). The body force cells are represented by the wire frame shaded
according to the strength of the axial component of the body force at each cell.
The axial velocities along the circumference of radius 0.7265 on this plane 0.323 radius upstream of the propeller are
shown in Figure 26. The velocities are shown at the moment the key blade is at the top position (0º). Since there are five
blades, the five velocity peaks are shown in the figure. The mean acceleration at this location is about 16% of the inflow
velocity, while the variation between blades is about 5.5% of the inflow velocity. The velocities predicted by the unsteady
Euler solver (WAKEFF-3U) and the potential flow solver (MPUF-3A) agree very well leaving only −0.4% to +0.7% error
in the effective wake velocity prediction.
Figure 25: DTMB 4842 propeller represented with body force cells; key blade at 0º.
5.4 Nominal Wake with Only the Third Harmonic
In this section, a nominal wake which only has a third cosine harmonic is considered with the zero-thickness DTMB
4148 propeller. Since the propeller has three blades and the inflow wake has only the third harmonic, the resulting
propeller loading should have dominant third harmonic components. The nominal wake has only the axial velocity which
can be described by the following equation.
(37)
The propeller working condition considered in this section is Js=0.8 and fully wetted flow. The fluid domain
covering −4≤x≤4 and 0.2≤r≤5 is discretized with 80x50x60 cells. Since the dimensionless time step size is determined to
be ∆t=0.004, 400 time steps correspond to one full revolution of the propeller. The applied artificial dissipation
coefficients are σ2=0 and σ4=1.
The convergence with revolution number is shown in Figure 27. The solution (u, v, w, p) at (x, r, θ)=(−0.323, 0.683,
the authoritative version for attribution.
0.0) when the key blade is at 0° (top) position of the ninth revolution is subtracted from that at the same location for the
same key blade position at the sixth, the seventh, and the eighth revolution. Since the resulting differences of the solutions
from those of the ninth revolution converge, the solution of the ninth revolution can be practically regarded to fully
developed; i.e. the solution repeats itself for consecutive revolutions.
2DEC Personal Workstation 500au is approximately equivalent to a 1 GHz Pentium PC.
lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
AN UNSTEADY THREE-DIMENSIONAL EULER SOLVER COUPLED WITH A CAVITATING PROPELLER ANALYSIS METHOD 632
Figure 26: Axial velocities predicted by the unsteady Euler solver and by the potential flow solver for the zero-thickness
DTMB 4842 propeller; xeff=−0.323.
Figure 27: The convergence with respect to the revolution for zero-thickness DTMB 4148 in the wake with only the third
harmonic.
Figure 28: The residual of the continuity equation and the propeller thrust and torque coefficients during the last (ninth)
revolution of zero-thickness DTMB 4148 in the wake with only the third harmonic.
Figure 28 shows the variation of the thrust and torque coefficients, KT and KQ, of the propeller during the ninth
revolution. As expected, the three propeller loading peaks, of which the amplitudes are about 50% of the mean loading,
are clearly observable. The curve labeled QM in the figure is the residual of the continuity equation, which is less than 10
−5 during the ninth revolution.
To observe the circumferential variation of the velocities, the velocities along the circumference at radius 0.7200 on
the plane x=−0.3 for the moment corresponding to the blade angle3, θ=0º, are plotted in Figure 29. Along the
circumference, the total velocity is varying between 0.59 and 1.08, while the propeller induced velocity varies from 0.8 to
1.4. Notice that the maximum of the propeller induced velocity occurs about 20º~30º after the blade. The effective wake
varies from 0.48 to 0.96 and this amplitude is smaller than that of the nominal wake which varies between 0.4 and 1.0. In
addition to these velocities, the steady effective wake computed by the steady Euler solver is plotted together. It is very
interesting that the unsteady effective wake at this moment (key blade at θ=0º) is very close to the steady effective wake.
In Figures 30 and 31, these velocities are compared for different moments corresponding to blade angles −24º and
−72º. Notice that the nominal wake and the steady effective wake (WAKEFF-3D) are, by definition, always the same. It
can be observed from the comparison of velocities at different mo
the authoritative version for attribution.
3A right-handed propeller is used in the present application. That is, the propeller rotates in the direction of decreasing θ.
lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
AN UNSTEADY THREE-DIMENSIONAL EULER SOLVER COUPLED WITH A CAVITATING PROPELLER ANALYSIS METHOD 633
ments that the unsteady effective wake does not vary very much with time, even though the total and the propeller
induced velocity do change. Also, it is found from this comparison that the unsteady effective wake is always very close
to the steady effective wake.
Figure 29: Circumferential variation of the velocities at r=0.72; Blade angle=0º.
Figure 30: Circumferential variation of the velocities at r=0.72; Blade angle=−24º.
Figure 31: Circumferential variation of the velocities at r=0.72; Blade angle=−72º.
Figure 32: Comparison of the mean and the third harmonic of the time averaged unsteady effective wake.
The time averaged unsteady effective wake can be calculated easily from the time series data of the unsteady
effective wake. In Figure 32, the circumferential mean and the third harmonic amplitude of the time averaged unsteady
effective wake are shown. The circumferential mean is about 0.03 greater than that of the nominal wake, while the third
harmonic amplitude is about 0.06 smaller than that of the nominal wake. In Figure 32, the mean and the third harmonic
amplitude obtained from the steady Euler solver (WAKEFF-3D) are shown together. The mean of the time averaged
the authoritative version for attribution.
unsteady effective wake is 0.721, at radius 0.72, which is very close to the mean of the steady effective wake with value
0.715, at the same radius. Also, the third harmonic amplitude of the time averaged unsteady effective wake is 0.241, while
that of the steady effective wake is 0.247. It can be observed again that the time averaged unsteady effective wake and the
steady effective wake are very close to each other for this case.
In conclusion, the unsteady effective wake is found to vary little with time for the case of DTMB 4148 propeller and
the nominal wake with only the third harmonic. As expected, the time averaged unsteady effective wake is also very
similar to the steady effective wake obtained by the steady Euler solver.
lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
AN UNSTEADY THREE-DIMENSIONAL EULER SOLVER COUPLED WITH A CAVITATING PROPELLER ANALYSIS METHOD 634
Figure 33: Predicted cavity shape, 60 circumferential cells.
Figure 34: Photographs of the cavity shape observed in CAPREX I.
5.5 Cavitation on Propeller DTMB 4148
Since it is well known that propeller cavitation is very sensitive to the non-uniformity of the inflow, the accuracy of
the predicted effective wake can be indirectly judged by the cavitation shape predicted for a given effective wake. The
cavitation can be predicted by the propeller potential flow solver once the effective wake is obtained through the
iterations between the Euler solver and the potential flow solver. It should be noted that in the evaluation of the effective
wake the propeller induced flow field has been determined by excluding the influence of the thickness and cavity source.
In this work, the cavitation shapes observed in CAPREX I (Mishima et al 1995) are considered. The effective wake
is predicted using the nominal wake measured in CAPREX II (Mishima et al 1995) and propeller DTMB 4148 working at
Js=0.9087 and the cavitation number based on the number of revolution, σn≡(pshaft−pvapor)/(0.5ρn2D2)=2.576 which is the
condition at which the photographs of cavitation were taken in CAPREX I. Note that the tunnel effect was included by
placing the circular tunnel wall boundary in the prediction of the effective wake. These effective wakes were predicted
using grids with 60 cells circumferentially. In Figure 33, the predicted cavity shapes are shown. The photographs taken in
the experiment CAPREX I are shown in Figure 34. The predicted cavity planforms are very similar to those shown in the
photographs.
Figure 35: Comparison of predicted axial velocity at r=0.68R with the measured one in CAPREX II.
5.6 Total Velocity in front of Propeller DTMB 4148
The velocity field predicted by the Euler solver can be directly compared with the velocity measured in CAPREX II
(Mishima et al 1995). In the experiment, the propeller DTMB 4148 was working at Js=1.05 in the non-axisymmetric
wake. The velocity measurement at 0.32 radius upstream of the propeller is selected to be compared with the predicted
velocity because this is the place where the effective wake is usually evaluated.
The velocity distribution along the circumference at radius 0.68 is shown in Figure 35, and that at radius 0.98 is
shown in Figure 36. The axial velocities predicted by the steady Euler solver (WAKEFF-3D), by the same steady solver
with the tunnel flow simulation (WAKEFF-3D, tunnel flow), and by the unsteady Euler solver (WAKEFF-3U) are
the authoritative version for attribution.
compared together. The velocity marked as “time average of WAKEFF-3U” is obtained by averaging the time history of
the velocity at a fixed point predicted by the unsteady Euler solver. The time averaging is exactly equivalent to the
method in which the measured data were processed during the experiment. In the unsteady Euler solution, the propeller
loading is obtained from the effective wake predicted by
lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
AN UNSTEADY THREE-DIMENSIONAL EULER SOLVER COUPLED WITH A CAVITATING PROPELLER ANALYSIS METHOD 635
the steady solver. This propeller loading includes the blade thickness effect because the purpose is to obtain the total
velocity and the effective wake is not evaluated in this computation.
Figure 36: Comparison of predicted axial velocity at r=0.98R with the measured one in CAPREX II.
Even though the predicted profiles cannot follow the sharp changes in the measured velocity, the global agreement
between the computed and the measured velocity is very good. The effect of the tunnel in this case is found to be very
small because the propeller loading is very light. Also, it can be observed from the figures that the time average of the
unsteady velocity field is very close to the steady velocity field.
6 CONCLUSIONS
Robust steady and unsteady Euler equation solvers, coupled with a previously developed propeller potential flow
solver, are developed for the prediction of the steady and unsteady effective wakes. To the author's knowledge, the
present unsteady Euler solver is the first method that ever computes the unsteady effective wake with sufficient space and
time resolution. Through extensively two-dimensional studies, it is shown that the surface distribution method, with
sufficient number of cells, can simulate the flow field very accurately. In addition, the camberline pressure model is
suggested for coarser grid applications and used in the present unsteady three-dimensional Euler solver. The blade
thickness effect is identified as the major source of error (in the order of 2%) in the prediction of the unsteady effective
wake. Omitting the velocity influence of the thickness source has been found to improve the accuracy of the evaluated
velocities, suppressing the error under 1% of the inflow velocity.
The present method is applied to a nominal wake inflow which contains only the third harmonic. It is found that the
unsteadiness in the predicted unsteady effective wake is small, and that the time average of the unsteady effective wake is
very close to the steady effective wake predicted by the steady Euler solver method. This finding leads to the conclusion
that the three-dimensional steady Euler solver can be a computationally efficient, yet accurate, practical tool for the
prediction of effective wakes. The cavity shape predicted using the effective wake obtained by the present method
showed good agreement with the photographs taken during the experiment. The circumferential distribution of the axial
velocity computed by the present method agreed also very well with the one measured upstream of the propeller.
As a future research, an improved body force model which includes the mass source to represent the blade thickness
can be considered. The inclusion of the mass source would also allow for a consistent representation of unsteady cavity
shapes within the Euler equation solver. Also, a way of linking the unsteady effective wake with the propeller potential
flow solver should be devised in the future. This is necessary to investigate if a truly unsteady way of feeding the
unsteady effective wake into the propeller potential solver could lead to the amplification of the unsteadiness.
ACKNOWLEDGMENT
Support for this research was provided by the following companies and research centers: American Bureau of
Shipping, Daewoo Heavy Industries Ltd., David Taylor Model Basin, El Pardo Model Basin, Hamburgische Schiffbau-
Versuchsanstalt (HSVA), Hyundai Maritime Research Institute, Kamewa AB, Michigan Wheel Corporation, Rolla SP
Propellers SA, VA Tech Escher Wyss (former Sulzer-Hydro GmbH), Ulstein Propellers AS, Volvo-Penta of the
Americas, and Wärtsilä Propulsion.
REFERENCES
ANDERSON, J. 1995 Computational Fluid Dynamics: The Basics with Applications. McGraw-Hill, Inc., New York.
BRESLIN, J. AND ANDERSEN, P. 1994 Hydrodynamics of Ship Propellers. Cambridge University Press.
CHEN, B., STERN, F. AND KIM, W. 1994 Computation of unsteady viscous marine propulsor blade and wake flow. Proceedings: Twentieth
Symposium on Naval Hydrodynamics. August.
the authoritative version for attribution.
lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
AN UNSTEADY THREE-DIMENSIONAL EULER SOLVER COUPLED WITH A CAVITATING PROPELLER ANALYSIS METHOD 636
CHOI, J.-K. 2000 Vortical inflow—propeller interaction using an unsteady three-dimensional Euler solver. Doctoral dissertation, Department of Civil
Engineering, The University of Texas at Austin. August.
CHOI, J.-K. AND KINNAS, S. 1998a A 3-D Euler solver and its application on the analysis of cavitating propellers. Proceedings: Proceedings of the
25th ATTC. American Towing Tank Conference, Iowa City, Iowa, September.
CHOI, J.-K. AND KINNAS, S. 1998b Numerical water tunnel in two and three dimensions. Journal of Ship Research, 42, 2, June, pp. 86–98.
CHOI, J.-K. AND KINNAS, S. 2001 Prediction of Non-axisymmetric Effective Wake by a 3-D Euler Solver. Journal of Ship Research (to appear).
CHORIN, A.J. 1967 A numerical method for solving incompressible viscous flow problems. Journal of Computational Physics, 2, pp. 12–26.
COURANT, R., FRIEDRICHS, K. AND LEWY, H. 1967 On the partial difference equations of mathematical physics. IBM Journal, Mar, pp. 215–234.
DAMLE, S., DANG, T. AND REDDY, D. 1997 Throughflow method for turbomachines applicable for all flow regimes. Journal of Turbomachinery,
119, Apr, pp. 256–262.
HUANG, T. AND COX, B. 1977 Interaction of afterbody boundary layer and propeller. Proceedings: Symposium on Hydrodynamics of Ship and
Offshore Propulsion System. March.
HUANG, T. AND GROVES, N. 1980 Effective wake: Theory and experiment. Proceedings: 13th Symposium on Naval Hydrodynamics. October.
HUANG, T., WANG, H., SANTELLI, N. AND GROVES, N. 1976 Propeller/stern boundary layer interaction on axisymmetric bodies: Theory and
experiment. Technical report. DTNSRDC 76–0113. DTNSRDC. December.
KERWIN, J., KEENAN, D., BLACK, S. AND DIGGS, J. 1994 A coupled viscous/potential flow design method for wake adapted multi-stage, ducted
propulsors using generalized geometry. Trans. SNAME, 102.
KERWIN, J., TAYLOR, T., BLACK, S. AND MCHUGH, G. 1997 A coupled lifting-surface analysis technique for marine propulsors in steady flow.
Proceedings: Propellers/Shafting '97 Symposium. Soc. Naval Arch. & Marine Engnrs., Virginia Beach, VA, September, 1–15 (Paper No. 20).
KINNAS, S., GRIFFIN, P., CHOI, J.-K. AND KOSAL, E. 1998 Automated design of propulsor blades for high-speed ocean vehicle applications. Trans.
SNAME, 106.
KINNAS, S., CHOI, J.-K., KOSAL, E., YOUNG, J. AND LEE, H. 1999 An integrated computational technique for the design of propellers with
specified constraints on cavitation extent and hull pressure fluctuations. Proceedings: CFD99—The International CFD Conference. June 5–7.
LARSSON, L. 1997 The state-of-the-art of CFD for ship hydrodynamics. Proceedings: The International CFD Conference. May 29–31, 1–17 (Paper
No. 0).
MISHIMA, S., KINNAS, S. AND EGNOR, D. 1995 The CAvitating PRopeller Experiment (CAPREX), Phases I & II. Technical report. Department of
Ocean Engineering, MIT. August.
NI, R.-H. 1982 A multiple-grid scheme for solving the Euler equations. AIAA Journal, 20, 11, Nov, pp. 1565–1571.
PATANKAR, S. 1980 Numerical Heat Transfer and Fluid Flow. Hemisphere Publishing, New York.
RHIE, C. AND CHOW, W. 1983 Numerical study of the turbulent flow past an airfoil with trailing edge separation. AIAA Journal, 21, 11, Nov, pp.
1525–1532.
SHIH, W.-Z. 1988 A combined Euler equation/surface panel solution to the shear interaction problem of an open or ducted propeller. Doctoral
dissertation, Department of Ocean Engineering, M.I.T. October.
STANIER, M.J. 1998 The application of RANS code to investigate propeller scale effects. Proceedings: 22nd Symposium on Naval Hydrodynamics.
August 9–14, 199–211.
STERN, F., KIM, H., PATEL, V. AND CHEN, H. 1988a Computation of viscous flow around propeller-shaft configurations. Journal of Ship Research,
32, 4, pp. 263–284.
STERN, F., KIM, H., PATEL, V. AND CHEN, H. 1988b A viscous-flow approach to the computation of propeller- hull interaction. Journal of Ship
Research, 32, 4, pp. 246–262.
STERN, F., KIM, H., ZHANG, D., KERWIN, J. AND JESSUP, S. 1994 Computation of viscous flow around propeller-body configurations: Series 60
cb =0.6 ship model. Journal of Ship Research, 38, 2, pp. 137–157.
TODA, Y., STERN, F., TANAKA, I. AND PATEL, V. 1990 Mean-flow measurements in the boundary layer and wake of a Series 60 CB=0.6 model
ship with and without propeller. Journal of Ship Research, 34, 4, pp. 225–252.
the authoritative version for attribution.
lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
AN UNSTEADY THREE-DIMENSIONAL EULER SOLVER COUPLED WITH A CAVITATING PROPELLER ANALYSIS METHOD 637
DISCUSSION
J.Kerwin
Massachusetts Institute of Technology, USA
This paper presents an extremely careful and systematic approach to the problem of hull-propeller interaction. I
agree with the authors that getting the effective wake right is particularly important in the case of predicting unsteady
propeller cavitation, since it is known that small changes in the inflow can result in large changes in cavity volume. The
authors' work is therefore particularly important since the weak link in the overall process of predicting both unsteady
cavitation and propeller-induced hull vibratory excitation is the uncertainty in the inflow.
As the authors point out in their introduction, we have been working for some time on developing the coupling
concept for both steady and unsteady propeller flows. In the case of steady flow, the global flow solver that we employed
was initially a steady, axisymmetric RANS code. While this procedure works well, we found that almost all of the
computing time was devoted to the RANS solution, with only a small part being associated with the vortex-lattice blade
solution. In addition, RANS gridding issues proved to be the major source of time delays in obtaining a solution.
We therefore decided to explore the feasibility of using MTFLOW, an Euler code developed by Professor Mark
Drela at MIT, in place of a RANS solver. This procedure is therefore closer to the steady-flow method presented by the
authors. Our experience so far has been that both computing times and user-preparation times have been substantially
reduced, while the accuracy of the solution is comparable.
Figure 1 shows a preliminary result obtained for a mixed-flow waterjet pump as part of an ongoing Naval Engineer's
thesis being carried out by C.J. Hanson. The circumferential mean streamlines are plotted, together with contours of swirl
velocity.
Of principal interest is the comparison of rotor thrust and torque predicted by the Euler code and the RANS code.
The results are quite close, in spite of the fact that this is a particularly complex flow case. The reduction in computing
time from 12 hours to 30 minutes is very substantial.
Figure 1. Comparison of results for a mixed-flow waterjet using a coupled Euler/Lifting surface code (MTFLOW/PBD) and a
coupled RANS/Lifting surface code.
I am particularly interested in the authors' comparison of their computed time-averaged effective wake and their
computed steady effective wake as shown in Section 5.4. We make the assumption, at the outset, that an axisymmetric
nominal wake remains axisymmetric. Thus, in the steady flow case, one only needs to couple the lifting-surface code with
an axisymmetric RANS or Euler solver. The authors' conclusion would seem to justify this assumption. Could the authors
elaborate on whether this assumption is valid for all practical cases, or whether there are situations where a full 3-D
unsteady Euler solution is required?
AUTHOR'S REPLY
The authors would like to thank Prof. Kerwin for taking time to review this paper and for making valuable comments.
It is very encouraging that the resulting thrust predicted by the coupled Euler/lifting surface code is quite close to the
result obtained by a coupled RANS/lifting surface code. The fact that the propeller thrust and torque predicted by the two
methods are quite close implies the closeness of the predicted effective wake inflows. When we decided to solve the Euler
equations, we neglected the effect of viscosity in the interaction between the vortical inflow and the propeller potential
flow, based on the following two assumptions. (a) The effect of viscosity is limited very near the wall boundary such as a
hub, and (b) the major player in the interaction is the inviscid vorticity dynamics. These assumptions seem to be valid
even for a very complex flow inside a mixed-flow water-jet according to Prof. Kerwin's numerical experiment.
the authoritative version for attribution.
lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
AN UNSTEADY THREE-DIMENSIONAL EULER SOLVER COUPLED WITH A CAVITATING PROPELLER ANALYSIS METHOD 638
The issue of axisymmetric analysis versus non-axisymmetric analysis of the effective inflow should be separated
from the discussion of the steady versus unsteady effective wake. As a final result of the current unsteady non-
axisymmetric Euler solver, we obtain the unsteady non-axisymmetric effective wake inflow, while the steady
non-axisymmetric effective wake, can be obtained from the iteration between the steady non-axisymmetric
Euler solver and the lifting surface code. Note that the body force used in this steady analysis is the time averaged
propeller force, which is varying in space. There is yet another prediction by using an axisymmetric Euler solver, in
which the velocity flow field is a function of only the radius, In the axisymmetric analysis, the body force is
obtained by circumferentially averaging the propeller force. One of the conclusions in the paper is that the time average of
is close to for the case tested in the current work.
In the course of our development of the axisymmetric and non-axisymmetric Euler codes, we were also interested in
the question if the circumferential average of would be close to In the several cases we have tested, we
have found that the circumferential average of is close to Therefore, we believe that the axisymmetric
analysis is also quite useful in predicting the mean performance of propulsors. Finally, we cannot yet generalize our
conclusion that the time average of is close to because we have very limited number of test cases.
One obstacle to reach a final conclusion on this issue is the fact that we cannot input the unsteady effective inflow into the
current lifting surface analysis due to the time invariant inflow assumption built in the method. In any case, the unsteady
Euler solver is still useful for inherently unsteady flow cases such as acceleration, deceleration, and maneuvering of a ship
the authoritative version for attribution.
