Twenty-Third Symposium on Naval Hydrodynamics (2001)

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

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. 

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, 

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 

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 

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. 

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 

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 

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. 

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. 

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). 

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 

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 627 with respect to the time it takes to travel a distance equal to the radius of the disk with a speed equal to the inflow velocity. Since the disturbance initiated by the sudden existence of the actuator disk travels with a speed faster than the dimensionless inflow velocity 1.0, it takes less than t=4 for this front of the disturbance to travel downstream and hit the downstream boundary. This is the reason why the changes between time steps are not decreasing as much until the dimensionless time reaches t≈3. Once the front of the disturbance hits the downstream boundary (at t≈3), the flow field is more or less fully developed, and thus the changes begin to decrease to reach a certain level limited by the machine accuracy. Figure 16: Convergence history of the actuator disk problem. Figure 17: Pressure and velocity evolution along the center line (curves are drawn at the dimensionless time interval of 0.5). Figure 18: Time evolution of the u-velocity field; t=1.0, t=2.0, t=3.0, t=4.0 (from top to bottom). In Figure 17, the time evolution of the axial velocity and the pressure is shown. Since the curves are drawn at the dimensionless time interval of 0.5, it can 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 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 

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<x<4 and 0.2<r<5. The dimensionless time step size for the unsteady solution is ∆t=0.004, and the artificial dissipation coefficients are σ2=0 and σ4=1. Figure 21: 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 one-bladed DTMB 4148 propeller. The axial velocities along the circumference at r=0.72 and axial location x=−0.323 at the instance the blade is at 0° position (top, or on y-axis) are shown in Figure 20. The velocity computed by the unsteady Euler solver (WAKEFF-3U) agrees very well with the velocity computed by the potential solver (MPUF-3A) at the angular positions away from the blade. The agreement between the two velocities becomes worse at the blade locations within ±50°, where flow deceleration and acceleration occur when the leading edge and the trailing edge of the blade, respectively, are aligned with the considered point. This is attributed to the fact that the camberline pressure model can not simulate accurately the effect of blade thickness on the flow near the blade, as was also demonstrated in section 4.1.2. The disagreement of the velocities results to an error in the effective wake ranging from −1% to 2%. In Figure 21, the three velocity fields at a location 0.323 radius upstream of the propeller are shown in order to observe the distribution of the total, the propeller induced, and the effective wake velocities on this plane. The velocities are plotted at the moment when the key blade is at the 0° (top) position. While the total and the propeller induced velocities agree very well globally, the two velocities near the blade differ by 1 to 2%. This difference results in the locally erroneous effective wake, as shown at the right of Figure 21. the authoritative version for attribution. The results shown up to now are computed for the propeller with the original thickness. To see the effect of blade thickness, a similar computation is performed with a zero-thickness propeller. The axial ve 

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.