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 509

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
SESSION 10
LIFTING-SURFACE FLOW: INVISCID METHODS

OCR for page 509

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

OCR for page 509

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
The Nonlinear Numerical Prediction of Unsteady Sheet Cavitation for Propellers of Extreme Geometry
N.E.Fine (Engineering Technology Center, USA)
S.A.Kinnas (Massachusetts Institute of Technology, USA)
Abstract
The unsteady flow around cavitating marine propellers is treated in nonlinear theory by employing a low-order potential based boundary element method. The solution is found in the time domain. The kinematic and dynamic boundary conditions, which are fully three dimensional, nonlinear and time-dependent, are satisfied on an approximate surface consisting of the propeller surface beneath the cavity and the portion of the blade wake surface overlapped by the cavity. An efficient and robust algorithm is developed to predict arbitrary cavity shapes, including so-called “mixed” cavity planforms in which the blade is partially cavitating at inner radii and supercavitating near the tip. In previous works, it has been shown that the present solution represents the first iteration of a completely nonlinear solution in which the exact boundary conditions are satisfied on the exact flow boundary. The main emphasis of the present work is to investigate the performance of the method for extreme propeller geometries, such as those with large amounts of skew, rake or pitch. The effect of the hub on the solution is investigated and shown to have more effect on heavily pitched propellers than on moderately pitched ones. The importance of including the crossflow terms in the dynamic boundary condition is also ivestigated for extreme propeller geometries.
1 Introduction
With recent increases in the demand for heavily loaded efficient marine propellers, the occurrence of cavitation has become less and less avoidable. As a result, it has become the task of the hydrodynamicist to predict the cavitation characteristics of a propeller at the design stage, with the expectation that analytical methods can be used to avoid excessive cavitation within an appropriate range of ship speed. Computational tools for the reliable prediction of propeller cavitation are therefore in high demand. Moreover, these tools must be applicable to propellers of extreme geometry (high skew, rake or twist) which have become very common in recent years.
The approach taken in this work is to treat propeller cavitation strictly as sheet cavitation. The travelling bubbles and bubble clouds which often are entrained in the wake of the sheet cavity, or precede the occurrence of sheet cavitation in the incipient stage, are not included in the present analysis. Viscosity will also be neglected in this work, and the inflow will be assumed to be free of vorticity. The advantage of this approach lies in its relative mathematical simplicity. The governing equation is Laplace's equation for the perturbation potential. The use of surface singularities and the now classic application of Green's 3rd identity are well suited for such a problem.
Despite the relative simplicity, however, the problem is far from trivial. The main difficulty arises from the need to determine the cavity free-surface on which the pressure is prescribed. The fact that a portion of the flow boundary is unknown makes the problem nonlinear. In addition, the inherent complexity of solving three dimensional flows around extreme geometries adds to the difficulty of the analysis. Furthermore, the cavity extent and volume variation in time resulting from nonaxisymmetries in the inflow is strongly nonlinear, making it impossible to solve in the frequency domain.
The main disadvantages of treating only sheet cavitation in potential flow result from the combined neglect of viscosity, tip and hub vortex cavitation, bubble cavitation, and cloud cavitation.

OCR for page 509

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
However, there are two additional rationales for taking the potential flow avenue. First, it is believed that the first-order contributor to dynamically varying blade loads and radiated pressures is the transient sheet cavity [28]. Second, many or all of the neglected phenomena may be included as refinements to the potential flow solution. For example, a viscous-inviscid boundary layer interaction could be used to model the effects of viscosity. As another example, a local tip solution, including a cavitating tip vortex, could be matched to the outer sheet cavity solution. The potential of such a model to capture most of the physics at reasonable computational expense is the motivation for obtaining an accurate and efficient potential flow solution which is able to treat a wide class of geometries including highly skewed, raked or pitched propellers. At this stage of development, it is also important to investigate the effects of the parameters which are included in the model. For example, the effect of the hub may be investigated to determine its importance to the solution for various geometries and operating conditions. Also, the importance of the crossflow terms in the dynamic boundary condition may be investigated. Finally, the method should be interrogated to determine if it predicts multiple solutions. These are the goals of the present work.
2 Formulation
The mathematical formulation first appeared in [15]. The formulation will be repeated here, but with more detail on the dynamic boundary condition in the wake (section 2.2.2). Details of the numerical implementation are provided in the next section.
Consider a cavitating propeller subject to a nonuniform inflow UW(xs,rs,θs), with the subscript 5 denoting the ship-fixed coordinate system in which the wake is defined. Uw is assumed to be the effective wake. Some details of the geometry are shown in Figure 1. The solution is found in the (x,y,z) coordinate system, which rotates with the propeller. The propeller is assumed (without loss of generality) to be right-handed and to rotate at a constant angular velocity ω. The inflow relative to the propeller is
(1)
Figure 1: The propeller and cavity geometry, coordinate systems, and nonuniform inflow. xs,ys,zs represent the ship-fixed coordinate system; x,y,z represent the propeller-fixed coordinate system. The point p is a point on the propeller surface, defined by the vector from the origin of the propeller-fixed system.
For the moment, assume that the propeller has a developed sheet cavity whose time-dependent surface is denoted by SC(t), as shown in Figure 1. The fluid is assumed to be inviscid and the resulting flow to be incompressible and irrotational. In this case, the time-dependent total flow velocity q(x,y,z,t), can be written in terms of the perturbation potential, (x, y,z,t), as follows:
q(x,y,z,t)=Uin(x,y,z,t)+(x,y,z,t). (2)
The goal is to determine the potential field (x,y,z,t), as well as the cavity surface SC(t)1. Once is known, the pressure distribution may be computed by numerically differentiating the potentials and applying Bernoulli's equation. The unsteady blade load distribution may then be determined by integrating the pressure. Knowing SC(t), the cavity volume history may be computed. In the following sections, the equations and boundary conditions necessary for the solution will be derived.
2.1 Green's Formula
The perturbation potential p(x,y,z ,t) at any point p which lies either on the wetted blade (or
1
The leading edge of the cavity will be assumed to be known, the rest of the surface Sc(t) is to be determined. See section 6 for related discussion.

OCR for page 509

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
Figure 2: Definition of the flow boundaries on which the boundary conditions should be satisfied.
hub) surface2, SW S(t), or on the cavity surface, SC(t), (both surfaces are shown in Figure 2) must satisfy Green's third identity:
(3)
The subscript q corresponds to the variable point in the integrations. is the unit vector normal to the surface of the propeller, the cavity, or the wake. The unit normal points into the fluid (on the wake surface, is oriented such that it points towards the same horizon as the normal on the suction side of the blade). Δw(r,θ,t) is the potential jump across the wake sheet, SW(t), and G(p; q)=1/R(p; q) is Green's function, where R(p; q) is the distance from the field point p to the variable point q.
Equation (3) expresses the perturbation potential on the surface formed by the union of the cavity and blade surfaces, SW S (t) ∪ SC(t), as a superposition of the potentials induced by a continuous source distribution, G, and a continuous dipole distribution, , on SWS(t)∪SC(t), and a continuous dipole distribution on the trailing wake surface, SW(t). The application of Green's 3rd identity to problems in potential flow is classic [18, 22, 26], while applications to propeller and rotor flows are more recent [21, 23, 7, 8, 10, 16].
Figure 3: The approximate cavity surface on which the boundary conditions are applied.
Ultimately, the exact nonlinear potential flow solution will be found when the kinematic and dynamic boundary conditions (which may be applied simultaneously with Green's formula (3)) are satisfied on the exact flow boundary. However, we face the usual problem that the position of the cavity surface is unknown. As a first iteration towards the fully nonlinear solution, we apply the boundary conditions on an approximate flow boundary. If the blade is experiencing only partial cavitation, then the approximate boundary coincides with the blade surface (including the part of the blade beneath the cavity). This surface will be referred to as SCB. On the other hand, if the blade is super-cavitating, the upper and lower sides of the super-cavity downstream of the blade trailing edge are collapsed to a single surface and the two sides of the collapsed cavity surface are taken to coincide with the two sides of the zero-thickness trailing wake sheet (see discussion below). This surface will be referred to as SCW(t). Thus the approximate flow boundary consists of the blade surface, SCB, and the portion of the trailing wake sheet which is overlapped by the cavity, SCW(t).
A sketch of the approximate boundary is shown in Figure 3. The geometry of the wake, SW(t) will from here on be assumed to be invariant in time and identical to the steady-flow relaxed wake
2
The wetted blade surface is that part of the blade which is not cavitating.

OCR for page 509

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
corresponding to the circumferentially averaged inflow [6]. The approximate flow boundary therefore coincides with that of the fully wetted solution, as described by Hsin [8].
A justification for making this approximation as well as a measure of its effect on the cavity solution is given in [3] and [15].
As a result of the linearization of the part of the supercavity downstream of the blade trailing edge, the surface SCW(t) may be considered to be a force-free surface. The pressure across the collapsed cavity surface must therefore be continuous and equal to the cavity pressure pc. The pressure across the blade wake surface must also be continuous. Therefore, it is possible to consider these two surfaces to coincide, with the condition of pressure continuity being
p+=p−=pc on SCW(t)
p+=p−=p on SW. (4)
Satisfying the boundary conditions on the approximate boundary described above may be viewed as the first iteration towards a fully nonlinear solution. In the fully nonlinear solution, subsequent iterations are found by satisfying the dynamic boundary condition on an updated cavity surface, where the kinematic boundary condition is used to update the surface, as described in section 2.3. The solution is then considered to be converged when the cavity surface does not change (to within a tolerance) between two consecutive iterations. However, it was discovered that, using the present potential based panel method, the first iteration solution (where Green's formula is satisfied on the approximate boundary) is extremely close to the converged nonlinear solution for a wide range of operating conditions [13,14,3]. As a result, it is deemed unnecessary to go beyond the first iteration towards the fully nonlinear solution. In view of the high computational cost involved in regridding and re-computing influence coefficients, the importance of this conclusion cannot be overstated.
Considering the approximate boundary, the first integral on the right-hand-side of equation (3) may be decomposed into integrals over the blade surface and the portion of the wake which is overlapped by the cavity (denoted SCB and SCW, respectively, in Figure 3). The exact form of Green's formula depends on the location of the field point, which will either be on SCB or on SCW. Each case will be considered separately in the following.
2.1.1 Field Point on SCB
If the field point is on the blade surface SCB, the local contribution is extracted from the first integral in equation (3) and Green's formula becomes
(5)
where the superscripts + and − correspond to the upper and lower sides of the wake surface, respectively. Note that the wake surface has been divided into SCW and SW, as shown in Figure 3.
The velocity normal to SCW is discontinuous across the surface. The jump in defines a source distribution, of density qw(t), which represents the cavity thickness:
(6)
The potential is also discontinuous across SW, where the jump Δw (r,θ,t) is related to the local circulation history.
Inserting (6) into (5) yields:
(7)
2.1.2 Field Point on SCW
Now consider the case where the field point is on the upper side of the collapsed cavity surface in the wake, SCW(t). Extracting the local contribution from the integral over SW on the right hand side of (3) and using the expression +(t)+−(t)=

OCR for page 509

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
2+(t)−Δw(r,θ,t) yields an expression for the potential on the upper wake surface +:
(8)
A detailed derivation of (8) may be found in [3]. Equations (7) and (8) define the potential (t) on the blade surface beneath the cavity, SCB, and the potential +(t) on the wake sheet, SCW(t), in terms of the following distributions of singularities:
continuous source and normal dipole distributions on SCB of strength (t) and (t), respectively
a source distribution on SCW(t) of strength qw(t)
a normal dipole distribution on the entire wake sheet SCW(t) ∪ SW of strength Δw(r,θ,t).
On the wake sheet SCW(t) ∪ SW the dipole strength Δw(r,θ,t) is convected along the assumed wake surface with angular speed ω, in order to ensure that the pressure jump in the wake is equal to zero:
(9)
where r,θ are the cylindrical coordinates of the wake surface and θT(r) is the θ coordinate of the blade trailing edge at radius r. is the unsteady fully wetted flow potential jump in the wake. is known since the unsteady fully wetted flow solution is found before the unsteady cavity solution.
Note that the wake convection given by (9) is identical to the fully wetted convection [17]. This is allowed by the linearization of the supercavity in the wake, where the two force-free surfaces are considered as one.
Everywhere on SCB and SCW(t), either the source distribution is known (a Neumann boundary condition) or the dipole distribution is known (a Dirichlet boundary condition). Green's formulae (7) and (8) may then be discretized and rewritten as a linear system of equations by employing the boundary conditions and shifting all of the known quantities to the right-hand-side.
2.2 Dynamic Boundary Condition
The dynamic boundary condition (DBC) requires that the pressure everywhere inside and on the cavity be constant and equal to the known cavity pressure, pc. Bernoulli's equation with respect to the propeller fixed system is
(10)
where ρ is the fluid density and r is the distance from the axis of rotation. Here qt is the total velocity on the cavity surface. po is the pressure far upstream on the shaft axis; g is the gravitational constant and ys is the vertical distance from the horizontal plane through the axis of rotation, as shown in Figure 1. ys is defined as negative in the direction of gravity.
After some manipulation, and using the definition of the cavitation number:
(11)
where and D are the propeller revolutions per unit time and diameter, respectively, we find that the magnitude of the total cavity velocity satisfies
(12)
The function f(s) corresponds to a pressure recovery law at the trailing edge of the cavity along the arc s on the surface of each span wise blade section. The pressure recovery is given by an algebraic expression over a specified portion of the cavity near its trailing edge. This termination model is described in detail in [3].

OCR for page 509

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
Since the cavity boundary consists of two parts —one coinciding with the cavitating portion of the blade surface and the other with the cavitating portion of the wake surface—further derivation of the dynamic boundary condition will be considered on each surface separately.
2.2.1 DBC on the Blade Cavity
In addition to the expression (12), the cavity velocity qt may also be expressed in terms of the directional derivatives of the perturbation potential and the components of the inflow along the same curvilinear coordinates [10]. The coordinate system3 on the cavity surface consists of s (chordwise) and ν (spanwise), as shown in Figure 3:
(13)
with and being the unit vectors corresponding to the coordinates s and ν, respectively, and with being the unit normal vector to the assumed cavity. Us, Uν, and Un are the s, ν and n components of the relative inflow, Uin.
If s, ν and n were located on the correct cavity surface, then the normal velocity, , would vanish. However, this is not the case since the cavity curvilinear coordinates are approximated with those on the propeller surface. Nevertheless, in applying the dynamic boundary condition, the normal velocity is assumed to be vanishingly small. In the fully nonlinear scheme, the normal velocity vanishes as the solution converges. As shown in [3] it was found that leaving the normal velocity term out of the dynamic boundary condition has only a small effect on the solution.
Equations (12) and (13) may then be combined to form an equation which is quadratic in the unknown chordwise perturbation velocity, . The two solutions to the quadratic equation represent potential gradients which are positive and negative in the direction of increasing s. The positive root is chosen to ensure that the cavity velocity points in the direction of increasing s. We can now express in terms of the cavitation number, the inflow velocity, and the unknown derivatives and :
(14)
with θ being the angle between s and ν, as shown in Figure 3. Here |qt|2 is defined as in equation (12). Equation (14) is integrated once to form a Dirichlet boundary condition on :
(15)
The integral on the right-hand-side of (15) is determined by trapezoidal quadrature. Since (15) defines the strength of the dipole distribution on the cavity, it may be directly substituted into Green's formula 7.
According to the dynamic boundary condition (15), depends on both its span wise and time derivatives. These terms will be treated as knowns and will be updated in a time-stepping scheme which will be discussed later. The influence of the crossflow term was studied first for the case of partially cavitating 3-D hydrofoils and it was found that the global dependence of the solution on the crossflow term was small. This result will be shown in section 4.2, where it will also be demonstrated for the propeller solution. For a discussion of the convergence of the time-derivative with iterations and its effect on the solution, see [3] or [15].
2.2.2 DBC on the Wake Cavity
The dynamic boundary condition on the cavitat— ing portion of the wake, SCW, may also be written as a Dirichlet condition on +. In this case, consider the orthogonal system (s,u,n), shown in Figure 4. Assuming that s is the direction of the mean velocity, the total velocity on the upper side of the wake sheet may be written
3
In general non-orthogonal.

OCR for page 509

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
Figure 4: Velocity diagram on the wake surface.
The normal velocity will be omitted from the dynamic boundary condition, as it was from (13), with the same justification.
Applying Bernoulli's equation, which is used to define the total velocity on the cavity qt, we have
(16)
The dynamic boundary condition on SCW may thus be written
. (17)
Equation (17) may be integrated once to form a Dirichlet boundary condition on the potential on the upper wake surface, +:
(18)
The integral in (18) is computed by trapezoidal quadrature. Equation (18) defines the potential + on the upper side of the wake and may be directly substituted into Green's formula (8).
2.3 Kinematic Boundary Condition
Since the dynamic boundary condition is applied on the portion of the boundary which is encompassed by the cavity, the other boundary condition, namely the kinematic condition, may be used to determine the position of the actual cavity surface once the singularity strengths have been determined. In this section, the most useful form of the kinematic boundary condition (KBC) will
Figure 5: Definition of the cavity surface defined with respect to the blade surface.
be derived. As in the previous section, the cavity boundary will be divided into two zones which will be considered separately.
2.3.1 KBC on the Blade Cavity
The kinematic boundary condition on the cavity is the requirement that the substantial derivative of the cavity surface vanishes:
(19)
where n is the coordinate normal to the blade surface (with unit vector ) and h(s,ν,t) is the thickness of the cavity normal to the blade at the point (s,ν) at time t. Expressing the gradient in terms of the local directional derivatives
(20)
performing the dot product with qt (as defined in (13)) and substituting the result in (19) yields the following partial differential equation for the cavity thickness:
(21)
where
.

OCR for page 509

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
Figure 6: Definition of the cavity surface defined with respect to the wake surface.
2.3.2 KBC on the Wake Cavity
The kinematic boundary condition on the cavity surface in the wake may be derived in a similar fashion, except that now both surfaces of the supercavity must be considered
(22)
where g±(s,u,t) defines the upper and lower cavity surfaces, as shown in Figure 6. Note that (s,u,n) is an orthogonal system. V+ and V− are the total velocities on the upper and lower sides of the wake (also shown in Figure 6), respectively, and may be written
(23)
The upper and lower cavity surfaces, g(s,u,t)±, may be written
where C is the cavity camber in the wake and hw is the cavity thickness. The quantities g, C and hw are all taken along the normal to the wake surface. These quantities are also shown in Figure 7.
Expanding equations (22) we find that
Figure 7: Definition of the cavity camber and height for a supercavitating section of the propeller blade.
(24)
Taking the difference between the two equations in (24) and assuming that ŝ coincides with the direction of the mean velocity so that
(25)
yields
(26)
Here we have used the definition of the wake source strength (6) and the following equalities:
(27)
which readily follow from the assumption (25) and the fact that the free vorticity must follow the mean velocity vector.
To be consistent with the dynamic boundary condition, we assume that the span wise crossflow velocity is small. Thus, applying equation (16), and the assumption of small span wise slope of the camber C(s,u,t), the kinematic boundary condition (26) reduces to
(28)
Note that the cavity height on the blade and in the wake, both shown in Figure 7, are defined differently and so are given separate symbols.
The position of the cavity surface over the blade surface is determined by adding the cavity thickness h normal to the blade surface at the midspan

OCR for page 509

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
of the panel boundaries. In the wake, the cavity camber C(s,u,t) must first be determined. An expression for C(s,u,t) may be found by adding the two equations (24) and dividing through by two:
(29)
where Un is the inflow velocity normal to the wake sheet. Equation (29) is numerically integrated to determine the camber surface in the wake. At the trailing edge of the blade, the continuity of camber and thickness is imposed:
(30)
Here, the superscripts + and − denote just upstream and just downstream of the trailing edge, respectively. is the value of the camber just upstream of the trailing edge. It is determined by adding to the trailing edge along the blade normal (see Figure 7). The quantity is determined by interpolating the upper cavity surface over the blade at the trailing edge and computing its normal offset from the wake sheet. The upper and lower surfaces of the cavity in the wake are then determined by adding and subtracting half of the cavity thickness hw from the camber surface. This defines the cavity surface at the midspans of all the spanwise strips. The surface of the cavity at the strip boundaries are determined by interpolation and extrapolation.
2.3.3 KBC on the Wetted Blade
The kinematic boundary condition on the wetted portion of the blade, SWS, defines the source strength there in terms of the known inflow velocity:
(31)
where xq,yq,zq are the coordinates of the point q with respect to the propeller fixed system. As in the case of fully-wetted flow, this boundary condition may be directly substituted in Green 's formula.
2.4 The Kutta Condition
The Kutta condition requires that the fluid velocity be finite at the blade trailing edge. It was found necessary, in the fully wetted steady flow propeller solution, to satisfy a nonlinear Kutta condition which ensures pressure equality between the suction and pressure sides at the trailing edge [20]. The so-called pressure Kutta condition was later refined by Kinnas and Hsin [16], who developed an efficient iterative solution. In the present work, an extension of Morino's steady Kutta condition is applied. Development of a pressure Kutta condition for the cavitating propeller is left for the future.
The value of the dipole strength, ΔT(r,t), at the blade trailing edge at time t is
(32)
where and are the values of the potential at the blade trailing edge at radius r on the suction side and the pressure side, respectively. The potential jump there is also equal to the circulation Г at time t around the blade section at radius r. This condition is equivalent to requiring the shed vorticity from the blade trailing edge to be proportional to the time rate of change of the circulation around the blade (Kelvin's law). An extension of Morino's Kutta condition for steady flow [24], in which the potential jump at the trailing edge of the blade is simply replaced by the potential jump at the nearest control points, is applied.
3 Implementation
The objective of the numerical analysis is to invert equations (7) and (8) subject to the kinematic boundary condition (31), the dynamic boundary conditions (15) and (17), and the Kutta condition. These equations, however, assume that the cavity extent is known. Since it is not, an iterative solution will be employed. This will be described in detail later in this section. First, we will describe how Green's formulae and the boundary conditions are satisfied for a given guess of the cavity extent. To accomplish this, we first tag one blade with the label “key blade”. The solution at a given time step will be obtained only for the key blade, with the influence of the other blades corresponding to an earlier solution of the key blade. The key blade surface is discretized into N chordwise and M spanwise quadrilateral panels with the corners lying on the blade surface SCB and with the control points located at the panel centroids. An example of a discretized blade is shown in Figure 8. The source and normal dipole distributions on each panel are approximated with constant strength distributions.

OCR for page 509

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
conditions are given later.
For 3-D flows the disturbance potential at infinity approaches zero. Consequently Eq. (9) is simplified to i=0.
In this case, the 3-D doublet strengths are obtained from the equations
(11)
Because we use uniform doublets in the 3-D case, Eq. (11) represents a unique system of N equations and N unknowns.
In both the 2-D (Eq. (10)) and the 3-D (Eq. (11)) cases, the vorticity shed at the trailing edge is given by the Kutta condition. This discussion follows.
Kutta Condition
The Kutta-Joukowski condition implicitly accounts for viscous effects otherwise neglected in potential-flow theory. In the present numerical method this is enforced by means of a pressure condition at the trailing edge. Accordingly, the pressures on the upper and the lower surfaces are required to be equal. For 3-D panels, upper and lower surfaces are referred to each specific strip of panels in the chordwise direction. This condition does not restrain the trailing edge from having a variation of pressure in the time domain or spanwise along the trailing edge.
Using the unsteady form of the Bernoulli equation (Eq. (5)), one obtains
(12)
Since the impermeability boundary condition is applied on δΩ, the flow on the upper and lower surfaces is tangent (superscript s) to both surfaces.
For two-dimensional problems, the scalar relations,
(13)
allow Eq. (12) to be recast into the form,
(14)
where and Г is the circulation of the airfoil defined as positive in a clockwise direction. The total blade circulation Г is linearly dependent on the shedding vortex strength γw as explained in the next section. Equation (14) is nonlinear in γ because VT is not a constant. In the present study, Eq. (14) is solved iteratively to account for this nonlinearity. Both Yao et al. (12) and Kim and Mook (13) ignored this nonlinearity and applied a more restricted Kutta condition at the trailing edge by requiring γ1=γN+1=0 and placing a wake vortex of unknown strength there. As shown later, the nonlinear approach gives improved solutions.
For the present 2-D system, Eq. (14) provides one of the two additional equations needed for solving Eq. (10). A second additional equation is needed to obtain a unique γ-distribution for a lifting body. This condition is provided by requiring the tangential velocity gradients along the trailing edge on both the upper and lower surfaces to be identical,
(15)
Requiring the equality of velocity gradients on both surfaces further ensures the smooth merger of the two flows. When second-order backward differencing is used, Eq. (15) is transformed to
(16)
where κu=SN−1/SN, κℓ=S2/S1 and κ=SN−1(1+κu)/S2(1+κℓ). Hence Eqs. (10), (14) and (16) form a determinate system of nonlinear equations for determining γ for a 2-D lifting configuration. They are solved using an iterative scheme.
For 3-D flows Eq. (12) has to be enforced at every panel along the blade trailing edge in the spanwise direction. Consequently, the analog of Eq. (14) for 2-D flow becomes a set of nonlinear matrix equations in 3-D flows. This system again requires iterative numerical procedures. For the 3-D case, Eq. (12) becomes

OCR for page 509

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
(17)
The perturbed form of Eq. (17) can be cast as
(18)
is the mean velocity at the blade trailing edge, and Δp =pu−pℓ, is the pressure difference between the upper and lower surfaces. In order to drive Δp to zero we adjust the shed vortex strength. This provides a mechanism for vortex shedding in 3-D calculations.
Calculation of Shed Vorticity
Kelvin's theorem states that the total circulation of the fluid at any instant is conserved. This condition provides a mechanism for shedding vorticity into the wake. For 2-D calculations a uniformly distributed vortex segment with strength γw is generated at each time step in the wake adjacent to the foil trailing edge. The length Sw of the vortex segment is set equal to the distance the trailing edge moves between t-Δt and t. At a subsequent time step this uniform line vortex is replaced by a discrete concentrated vortex of equivalent strength located at the center of the segment. The generated trailing-edge vortex segment is related to the total foil circulation, Г, as
(19)
The model depicted in Eq. (19) yields numerically stable solutions that are insensitive to the time-step used. When the nonlinear pressure condition is selected for solving Eq. (14), representative calculations based upon the concentrated vortex model used by Kim and Mook (13) for modeling the trailing-edge vorticity generation were found to be dependent on the time step used
Following the shedding process at the trailing edge, the wake vortices are convected downstream and develop a vorticity field. In the 2-D form of the present numerical scheme, these vortices are tracked through the flow field using Lagrangian techniques. The convection of these discrete vortices is modeled by a predictor-corrector scheme, given in Ref. (2).
The shed vorticity distribution in 3-D flows is obtained by solving Eq. (18) iteratively by relating the nonzero Δp to Δγw by
(20)
where D and are both matrices representing the disturbed velocity potential and induced velocity due to ΔγW of unit strength. These matrices are only dependent on the blade geometries and are required to be calculated only once. Thus Eq. (18) becomes
(21)
The sub-iteration formula at each time step to update shed vorticity strength is then given by
(22)
Vortex Core Splitting and Merging
As Eq. (19) indicates, a vortex is generated at each time step. This mechanism of generating wake shed vorticity will produce tracing instability when the wake vortices start interacting with one another (14). A remedy for redistributing and splitting the vortex cores was proposed by Mook et. al. (14). In addition to the predictor-corrector tracing method mentioned above, the present vortex convection scheme also adopts Mook's splitting and merging techniques. When any pair of vortices is within 1/4 u o Δt, a merging process takes place to reduce the number of vortices. On the other hand, when the distance between two vortices becomes larger than 1.5 uo Δt, a linear splitting process is used. The total circulation around the entire wake and airfoil is, however, maintained exactly. This vortex core splitting procedure maintains reasonably uniform spacing between vortices.
Cascade Formulation
For a 2-D cascade of blades with pitch H, each blade generates circulation. If the cascade runs along the y-axis, there exists an upwash far upstream of the cascade and a downwash far downstream. In conjunction with a specified inflow onset condition, an extra term is needed to ensure a unique upstream inflow condition, i.e.

OCR for page 509

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
(23)
for a multi-blade-row cascade flow, where MJ is the total number of blade rows, Hmj is the mj-th cascade spacing or pitch, and Гmj is the blade circulation from the mj-th cascade.
A complete integration of the velocity potentials for the source and the vortex distributions on each segment was performed numerically for the cascade modification. The detailed formulation as well as the modification in Green function for a cascade can be found in Ref. (2).
Gauss-Seidel Solution Procedure
For a 2-D problem, the governing matrix equations are usually small in size, even for multi-row cascade and therefore a direct solver can be used to solve Eq. (10) efficiently. For a 3-D problem, however, the size of the matrix equations is at least an order of magnitude larger than that of the 2-D ones and often reaches a few thousands or more. Consequently, direct matrix solvers are too expensive to execute for matrix equations of this scope.
In the present paper a two-level multi-block Gauss Seidel procedure is developed for obtaining solutions of 3-D problems, First, the left-hand side of Eq. (11) is partitioned on the basis of the larger surface patches such as those on the blade and hub. Thus Eq. (11) becomes
(24)
where NPS(ℓ) and NPE(ℓ) are the starting and ending panel numbers on surface patch ℓ, RHS is the right-hand-side terms of Eq. (11), NS is the total number of surface patches, and is the latest value of γj. On the ℓf-th surface patch, γj(k+1) is updated by
(25)
where RHS*denotes the right-hand-side terms of Eq. (24), and i1 and i2 are indices of the two trailing-edge panels on the same strip of the blade. The advantage of Eq. (25) is that it retains the dominant terms in the left-hand-side matrix of Eqs. (24).
RESULTS AND DISCUSSIONS
The unsteady rotor-stator interaction includes both potential-flow and wake effects. The validation of the current calculation method for a complex rotor-stator system therefore involves two validations of fundamental flow phenomena. The first of these is to demonstrate that the shed vortices from a foil (or cascade) are accurately computed. The second is to demonstrate that the interaction between a vortex and a foil is properly computed. Both of these fundamental phenomena have been demonstrated in Ref. (7).
The capability for accurately modelling shed vorticity was demonstrated by comparing the predicted wake behind an oscillating airfoil with experimental measurements (15). Although the present inviscid approach can not predict the viscous (velocity defect) wake at low blade oscillation frequency, it predicts well the jet-like wake for high frequency oscillation. The vortex and blade interaction process was also demonstrated in Ref. (7) by the computation of the interaction between a free traveling vortex encountering a stationary airfoil (11). The results indicate that the current method is able to accurately predict the interaction forces acting on the foil due to a passing vortex.
A third phenomenon that occurs in rotor-stator interactions is wake vortex cutting. An Example of this was given by Chow and Huang ( 8) who presented an analytical solution based on conformal mapping for a close encounter between a Joukowski airfoil and a moving vortex. The flow phenomenon described by Chow assumes that a “neutral release ” position exists for the incoming vortex which divides the vortex trajectories into those either above or below the airfoil.

OCR for page 509

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
Figure 1. Predicted instability of a vortex interacting with a Joukowski foil
Representative calculations for vortices on either side of the neutral release point are presented in Figs. 1 and 2 along with corresponding solutions from Chow and Huang. Figure 1 compares the trajectories of the free vortex at representative time steps for both methods. Because Chow's analytical solution extends only a chord length downstream of the airfoil, the comparison is limited to this region. For perspective we also present a snapshot of the wake vortices at the last time step of our calculation. The curved trajectories wrapping around the airfoil and the large amplitude wake roll-up indicate a very strong interaction between the vortex and the airfoil and thus provide a more critical validation for the present calculation scheme than that in Ref. (7). The corresponding force comparisons for these calculations are given in Fig. 2. When the free vortex travels above the foil, it is first absorbed by the foil as the foil's lift increases. After the vortex passes the trailing edge, the foil repels the vortex and its lift is also dropped. The lift on the foil is eventually leveled off as the vortex travels far downstream. Similarly, when the vortex starts from a location above the foil and travels to the lower side of the foil 's surface, it is again absorbed by the foil. As the vortex gets too close to the foil surface, the surface vorticity of the foil starts repelling the free vortex. At this moment, the lift of the foil drops suddenly and the vortex is repelled backward and starts moving towards lower side of the foil. On the foil's lower side, the trajectory of the vortex is much closer to the foil surface than that on the upper-side motion. During this period, the lift of the foil increases as the vortex moves towards the trailing edge. When the vortex is in the vicinity of the trailing edge, it starts interacting with the shed wake vorticity originating from the foil trailing edge. The spike in the predicted lift distribution near the trailing edge is related to the this strong interaction and depends also upon the Kutta condition used. Nevertheless, the present predicted lift distributions compare favorably with Chow's solutions.
Figure 2. Comparisons of present predicted force on the Joukowski foil with Chow's solutions
An additional demonstration of the effectiveness of the present method, and one that is much more similar to a rotor-stator interaction, is the flapper and foil experiment recently carried out at MIT ( 9,18). In this experiment detailed measurements of the unsteady flowfield around a foil were measured using the setup shown in Fig. 3. In this experiment two small NACA 0025 flapping foils of 3-inch chord were configured

OCR for page 509

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
upstream of a larger stationary foil, referred to as the “foil” in the later text, of 18-inch chord. The flappers performed synchronized sinusoidal motions of 6-degree amplitude at a reduced frequency of 3.62. The flappers generate periodic wake fluctuations which impose an unsteady condition on the foil. Velocity and pressure measurements at a Reynolds number of 3.78×106 were taken on the foil and around the foil in a rectangular box downstream of the flappers and enclosing the foil, as shown by the dash line in the Fig. 3. The box measurements were intended to provide upstream and downstream boundary conditions for Navier-Stokes computations. Unsteady flow separation was observed on the suction side near the trailing edge of the foil.
Figure 3. Schematic of MIT Flapping Foil experiment
Figure 4. Six-row computational cascade for simulating MIT's flapper and foil experiment
To simulate this flow configuration, a series of image blades were added to the geometry shown in Fig. 3 to incorporate the effects of the tunnel walls. Because of the asymmetry of the foil, the complete system of flappers, foil and image blades is a six-blade-row cascade system as shown in Fig. 4. Since a thick boundary layer exists in the aft region of the foil, a displacement thickness obtained from boundary-layer calculations based on steady solutions to the six row cascade system was added to the foil surface. Computed steady foil surface pressures with and without the displaced thickness correction are shown in Fig. 5. Comparisons of steady pressure distributions for calculations with and without the flappers and its image cascades are nearly the same.
Figure 5. Predicted steady-state pressure distributions on MIT's foil
Figure 6. Predicted unsteady lifts on MIT's foil and flappers

OCR for page 509

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
The unsteady calculations were started at t=0 from the steady solution. Figure 6 shows the instantaneous lift distributions of the flappers and the foil plotted against nondimensional time at reduced frequencies of 3.62 and 7.2. The foil transient lasts about one period and the calculation extends for six periods with 50 time steps in each period. The corresponding time-averaged mean pressure distributions from the unsteady calculation (reduced frequency of 3.62) are shown in Fig. 7 along with the envelopes of the unsteady solution. These are compared with the measured means and measured local maxima and minima, indicated by I-bars. The calculated mean unsteady pressure agrees well with the measurements. Since the present calculation uses an open trailing-edge for the displaced foil, fluctuations in the predicted pressure distributions exist near the trailing edge, as seen in Fig. 7. Figures 8 and 9 present comparisons of the time history, amplitude and phase angle of Cp at several locations on the stationary foil. Overall, these comparisons in Fig. 8 show the unsteady effects are quantitatively reasonable, but the amplitude are significantly larger than the experimental results. Near the trailing edge the overprediction becomes worse. This is likely due to the aforementioned pressure fluctuation introduced by the open trailing edge (see comparison at x/c=0.972 on Fig. 8).
Figure 7. Predicted time-averaged unsteady pressure distributions and unsteady pressure envelopes on MIT's foil
This overprediction is further demonstrated in Fig. 9 by the amplitude of the first harmonic of the pressure distribution around the airfoil. Again the computed amplitude on the suction side (SS) are much larger than the experimental measurements. Close inspection of the measured time history on Fig. 8, however, shows that the experimental variation of the pressure is predominantly a double frequency variation that is twice the reported flapper frequency. For this reason, the calculation was repeated at the double frequency (7.2 reduced frequency) for the flappers. The history of this lift was previously shown in Fig. 6. As this figure shows the amplitude of the fluctuating lift on both flappers is about the same for both high and low frequencies, but that for the stationary foil the amplitude is reduced considerably at the higher frequency. This phenomenon is also evident in the time history and the fundamental harmonic (7.2 reduced frequency) of the pressure distribution (Figs. 8 and 9). Both figures show the predictions of unsteady pressures and its amplitudes. These results show that predictions at high frequency agree surprisingly well with measured values. This suggests a possible discrepancy in the measured frequency of the flapper motion. Corresponding comparisons of the phase distribution of the pressure fluctuation around the foil are also given in Fig. 9 for both reduced frequencies. The predictions show little effect of frequency on the suction side, whereas on the pressure side the predictions show a strong phase lead near the leading edge at the lower frequency. As usual in such comparisons, the phase angle predictions do not compare as well with experiments. This phase comparison thus provides no illumination of the above dilemma.
Figures 10 through 11 show predictions of the rotor-stator interaction for the UTRC single stage, low speed turbine (10). The turbine consists of a 22-blade stator with pitch H=0.854c and a 28-blade rotor with pitch H=0.813c. Both the stator and the rotor have round trailing edges. The gap between the two blade rows is 15% of the stator chord. The rotor operates at a reduced passing frequency of 2.8, based on the turbine diameter.
To reduce the computation effort the 22 stator blades versus 28 rotor blades were simulated at the blade number ratio of 21 stator blades versus 28 rotor blades. This enables us to use three representative stator blades and four rotor blades. Trajectories of the shed vorticity patterns at four different time steps during the start-up period of the stator-rotor interaction, are shown in Fig. 10. As the wakes develop from the rest condition they initially follow a straight line behind the trailing edges of both the stator

OCR for page 509

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
Figure 8. Time history of Cp at various locations of the foil

OCR for page 509

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
Figure 9. Comparisons of amplitude and phase angle of the first harmonic of Cp
Figure 10. Snapshots of predicted vortex patterns at each time step for UTRC's stator-rotor turbine.
Figure 11. Pressure predictions on UTRC's stator-rotor turbine

OCR for page 509

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
and the rotor (t=10). At later times the stator wake begins to roll up as it enters the moving rotor (t=100). This roll-up of the stator wakes continues, and as they are cut by the rotor they begin to approach a periodic state (t=300 and t=500). Throughout, the wakes from the rotor remain essentially straight.
The time-averaged pressure distribution from one period of the later stages of the calculation in Fig. 10 are given in Fig. 11 along with the mean experimental measurements for both the stator (Fig. 11b) and the rotor (Fig. 11d). For reference we also show passage-averaged pressure distributions from quasi-unsteady stator-rotor calculations (Figs. 11a and 11c). These quasi-unsteady calculations were performed as a series of steady state calculations with various relative positions for the rotor-stator blades. Wake vortex shedding was neglected.
The averaged pressure distributions from both the quasi-unsteady calculation and the fully unsteady calculations are in good agreement with experimental measurements. It suggests that for preliminary estimates, especially at an initial stage of design when a large number of comparison studies are needed in a short time period, quasi-unsteady computations should be appropriate for providing reasonable performance predictions. The envelopes of the maxima and minima for both the unsteady and the quasi-unsteady calculations are shown in Fig. 11 along with the corresponding experimental measurements. This comparison clearly shows that the fully unsteady computation which includes vortex shedding and a convection effect has a distinctive advantage over the quasi-unsteady result. First the time-dependent loading is spread evenly over most of the blade. This even spreading of the unsteady loading over the entire blade surfaces is confirmed by the experimental measurements. The quasi-unsteady calculations predict the maxima and minima, variations are restricted to locations near the trailing edge of the stator and the leading edge of the rotor, and the magnitudes of the fluctuations in these regions are very large, far exceeding the experimental measurements. The primary reason for this even spreading of the pressure fluctuations in the true unsteady calculation is the presence of the wake vortex field. The wake vortices are continuously cut by blades as they are convected downstream through the rotor passages. These rolled up vortices have a significant impact on the unsteady pressure distribution on both the rotor and the stator. Clearly, the computational structure of the wake rollup (see Fig. 10) would be impossible to describe in advance.
CONCLUSIONS
A time-accurate hybrid panel/Euler prediction method is presented for unsteady incompressible flows through a rotor-stator system. By tracing wake vorticity, this method is capable of predicting complex flow features such as wake cutting and wake wrapping in rotor-stator systems.
The time-averaged blade pressure distribution predicted by the present approach is in good agreement with experimental measurements. Current results also suggest that for rotor-stator systems quasi-unsteady approaches are appropriate for time-mean performance predictions.
For time-dependent solutions of multiple-blade-row cascade flows, results obtained from the present fully unsteady approach compare well with experimental measurements. Details of the flowfield, such as wake vortex shedding, wake-wrapping and wake-cutting, are simulated for rotor-stator systems from the flow transient to the fully developed periodic stage. For rotor-stator systems quasi-unsteady approaches overpredict the unsteady loading in the strong interacting area of the blade and considerably underpredict in the weak interacting area. By including the complex vortex convection in the flowfield, the present method predicts a much more evenly distributed unsteady pressure as the measurement indicates.
Calculations on the MIT's flapping foil simulation have suggested that the multiple blade row cascade calculation method has correctly modelled the unsteady water-tunnel flow. The unsteady lift on the flappers is independent of its vibrating frequency, but the amplitude (the mean is the same) of the lift on the stationary foil reduces as the flappers move faster. The predictions based on the higher frequency flapper motion, however, agrees better with the measurements in local pressure distributions.
Our results also suggest that the multi-row blade interactions in many rotor-stator systems are inviscid in nature and can be analyzed with more efficient and inexpensive approaches, such as the present method.
The present time-accurate approach can be applied to other unsteady flow problems such as blade passage noise predictions and unsteady cavity flow predictions in rotor-stator systems.
ACKNOWLEDGEMENTS
This work was supported by the Independent Research Program and also by the Submarine Block Program for Hydrodynamics and Hydroacoustics of Internal Flow administrated at the David Taylor Model

OCR for page 509

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
Basin. Partial computational results were obtained using the facilities sponsored by the NASA Ames Numerical Aerodynamic Simulation Program.
REFERENCES
1. Giles, M.B., “Calculation of Unsteady Wake/Rotor Interaction,” AIAA J. of Propulsion and Power, Vol. 4, No. 4, 1988, pp. 356–362.
2. Lee, Y.T., Bein, T.W., Feng, J., and Merkle, C.L., “Unsteady Rotor Dynamics in Cascade,” ASME J. of Turbomachinery, Vol. 115, 1993, pp. 85–93.
3. Kerwin, J.E., and Lee, C.S., “Prediction of Steady and Unsteady Marine Propeller Performance by Numerical Lifting-Surface Theory,” SNAME Transactions, Vol. 86, 1978, pp. 218–253.
4. Preuss, R.D., Suciu, E.M., and Morino, L., “Unsteady Potential Aerodynamics of Rotors with Applications to Horizontal-Axis Windmills,” AIAA Journal, Vol. 18, No. 4, 1980, pp. 385–393.
5. Idris, B.M., Maruo, H., and Ikehata, M., “Theoretical Analysis of Unsteady Characteristics of Marine Propeller in Ship's Wake,” J. of the Society of Naval Architects of Japan, Vol. 156, 1984, pp. 60–68.
6. Shoji, H., and Ohashi, H., “Lateral Fluid Forces on Whirling Centrifugal Impeller (1st Report: Theory,” Transactions of the ASME, Vol. 109, 1987, pp. 94–109.
7. Lee, Y.T., Feng, J., and Merkle, C.L., “Prediction of Vortex and Linear Cascade Interaction Noise,” ASME Paper 93-GT-314, 1993.
8. Chow, C.Y., and Huang, M.K., “Unsteady Flows About a Joukowski Airfoil in the Presence of Moving Vortices,” AIAA Paper No. 83–0129, 1983.
9. Kerwin, J., Keenan, D., Mazel, C., Horwich, E., and Knapp, M., “MIT/ONR Flapping Foil Experiment: Unsteady Phase,” MIT/ONR FFX Workshop, David Taylor Model Basin, March 1993.
10. Dring, R.P., Joslyn, H.D., Hardin, L.W., and Wagner, J.H., “Turbine Rotor-Stator Interaction,” ASME J. of Engineering for Power, Vol. 104, 1982, pp. 729–742.
11. Lee, D.J., and Smith, C.A., “Distortion of the Vortex Core During Blade/Vortex Interaction,” AIAA Paper 87–1243, 1987.
12. Yao, Z.X., Garcia-Fogeda, P., Liu, D.D., and Shen, G., “Vortex/Wake Flow Studies For Airfoils in Unsteady Motions,” AIAA Paper 89–2225, 1989.
13. Kim, M.J., and Mook, D.T., “Application of Continuous Vorticity Panels to General Unsteady Incompressible Two-Dimensional Lifting Flows,” J. of Aircraft, Vol.23, No.6, 1986, pp. 464–471.
14. Mook, D.T., Roy, S., Choksi, G., and Dong, B., “Numerical Simulation of the Unsteady Wake Behind an Airfoil,” J. of Aircraft, Vol. 26, No. 6, 1989, pp. 509–514.
15. Koochesfahani, M.M., “Vortical Patterns in the Wake of an Oscillating Airfoil,” AIAA J., Vol. 27, No. 9, 1989, pp. 1200–1205.
16. Fuhs, D., “Comparisons of Calculations and Measurements,” MIT/ONR Flapping Foil Workshop, Washington D.C., March 29–30, 1993.

OCR for page 509

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
DISCUSSION
by Dr. S.Kinnas, MIT
The authors presented an interesting application of inviscid flow methods for practical design. It would be very interesting also if they show convergence tests of the involved methods, in particular, with number of panels and time step size for different values of reduced frequency of the incoming gust.
Author's Reply
The numerical convergence for steady flow calculations has been reported in reference (2) of the main text. Numerical errors, defined as the root-mean-square values of the differences between the calculated Cp and the exact Cp distributions, were presented based on different panel numbers for both circular cylinder and NACA 0010 airfoil. It was concluded that the present approach generates results with an accuracy better than the second order methods and the panel number in the order of 100 for airfoil geometries is adequate to drive the error down to 10 −5.
For unsteady flows, the panel size is also governed by two other factors. First, the size of the panels near the foil trailing edge should not be larger than the distance between two successive vortices shed from the trailing edge.
That is
Δs=Δt•Uo (A1)
where Δt=2π/ωnt and nt is the number of time steps in one period. And Eq. (A1), in terms of reduced frequency f=ωc/2uo, becomes
(A2)
Second, when downstream foils cutting through a vortical wake shed from upstream foils, the size of the panels adjacent to the wake vortices is limited to the distance traveled by the vortices in one time step during the close encounter between the vortex and the foil. This condition requires a similar constrain as shown in Eq. (A2) for panels in the front one-third of the foil.
For a typical calculation with a reduced frequency of 10, a normalized unity chord length and freestream velocity, and 25 time steps in one calculation period, c/Δs is about 80. On the average, we use 25 to 50 time steps and 160 to 300 panels to represent periodic foil motion in unsteady flows.