Abstract:
A method for calculating parameters about an axisymmetric body in a cavity is provided. The user provides data describing the body, a cavity estimate, and convergence tolerances. Boundary element panels are distributed along the body and the estimated cavity. Matrices are initialized for each panel using disturbance potentials and boundary values. Disturbance potential matrices are formulated for each panel using disturbance potential equations and boundary conditions. The initialized matrices and the formulated matrices are solved for each boundary panel to obtain panel sources, dipoles and cavitation numbers. Forces and velocities are computed giving velocity and drag components. The cavity shape is updated by moving each panel in accordance with the calculated values. The method then tests for convergence against a tolerance, and iterates until convergence is achieved. Upon completion, parameters of interest and the cavity shape are provided. This invention also allows determiniation of cavity shape for a cavitation number.

Description:
STATEMENT OF GOVERNMENT INTEREST 
   The invention described herein may be manufactured and used by or for the Government of the United States of America for governmental purposes without the payment of any royalties thereon or therefore. 
   CROSS-REFERENCE TO RELATED PATENT APPLICATIONS 
   Not applicable. 
   BACKGROUND OF THE INVENTION 
   (1) Field of the Invention 
   The present invention relates to computer model of hydrodynamic flows and more particularly, relates to modeling partially cavitating flows over a supercavitating axisymmetric body. 
   (2) Description of the Prior Art 
   Modeling of boundary flows about objects subject to laminar and turbulent flows is well known in the art. High speed underwater vehicles, however, cause cavitation of the surrounding fluid. Cavitation reduces pressure in the fluid below its vapor pressure causing the fluid to vaporize, allowing the undersea vehicle to travel with lower friction when the vehicle is completely surrounded by the cavity. 
   Partial cavitation is an unsteady phenomenon that occurs when part of the supercavitating vehicle is traveling in the cavity. Specifically, this phenomenon occurs during launch of the vehicle. A steady, partial cavitation allows development of vehicle designs which take advantage of drag reduction through cavitation. It may also be possible to take advantage of drag reduction with partial cavitation by properly directing the re-entrant jet that forms in the cavity closure region. Partial cavitation often occurs during maneuvering of the supercavitating vehicle. 
   A slender body theory has been developed to solve axisymmetric supercavitating flows. Using the slender body method, sources are defined along the body-cavity axis and control points along the body-cavity surface. A nonlinear differential equation is formed by imposing dynamic boundary conditions on the cavity. A conical cavity closure is assumed in order to solve the developed nonlinear differential equation. 
   A non-linear boundary element method for determining a cavity shape has been developed. Source and dipole strengths along the body-cavity surface are determined using kinematic boundary conditions on the wetted body surface and dynamic boundary conditions on the assumed cavity shape. The kinematic boundary condition is then used to update the cavity shape. The process is then iterated to solve for the unknown cavity shape. 
   Two numerical hydrodynamics models have been developed by the Naval Undersea Warfare Center for axisymmetric super cavitating high speed bodies. These models are the slender body theory (SBT) model and the boundary element (BE) model. Both of these models have been proven to predict cavity shape and parameters with good accuracy. 
   These models, however, do not account for the transition case when the vehicle is subjected to only partial cavitation. 
   In the SBT model, total drag is predicted by adding the pressure drag obtained from the model solution and the viscous drag obtained by applying the Thwaites and Falkner-Skan approximations along the wetted portions of the cavitator. This method is extended to subsonic compressible flows using the compressible Green&#39;s function. In the BE model, sources and dipoles are defined on the body-cavity shape and are solved using Green&#39;s formula. This yields a Fredholm integral equation of the second kind which gives the supercavitating cavity shape. 
   Partial cavitation modeling has been done by Uhlman, J. S. (1987), The Surface Singularity Method Applied to Partially Cavitating Hydrofoils, Journal of Ship Research, Vol. 31, No. 2, pp. 107-24; Uhlman, J. S. (1989), The Surface Singularity or Boundary Integral Method Applied to Supercavitating Hydrofoils, Journal of Ship Research, Vol. 33, No. 1, pp. 16-20; Kinnas, S. A., and Fine, N. E. (1990), Non-Linear Analysis of the Flow Around Partially and Super-Cavitating Hydrofoils by a Potential Based Panel Method, Proceedings of the IABEM-90 Symposium, International Association for Boundary Element Methods, Rome, Italy, and Kinnas, S. A., and Fine, N. E. (1993), A Numerical Nonlinear Analysis of the Flow Around Two- and Three-Dimensional Partially Cavitating Hydrofoils, Journal of Fluid Mechanics, Vol. 254. However, these methods are explicitly adapted for hydrofoils, and the theories presented therein are not readily adapted to supercavitating vehicles. 
   SUMMARY OF THE INVENTION 
   One object of the present invention is a method for modeling partial cavitation. 
   Another object is that such method model partial cavitation about an axisymmetric vehicle. 
   Accordingly, the present invention provides a method for calculating cavity shape for partial cavities about an axisymmetric body having a cavitator located at the foremost end. The method includes receiving system parameter data including geometric data describing the axisymmetric body, a cavity length, and a convergence tolerance. Boundary element panels are distributed along the body-cavity surface and matrices are initialized for each boundary element panel using the unit dipole, unit source functions and known boundary values. Disturbance potential matrices are formulated for each boundary element panel using disturbance potentials, normal derivatives of disturbance potentials, and no net flux boundary conditions. The initialized matrices and the formulated matrices are solved for each boundary panel to obtain unknown disturbance potentials along the wetted body-cavity surfaces, and normal derivatives of disturbance potentials along the cavity surface. The cavity position is then updated by moving each panel to satisfy the kinematic boundary condition, no flux across the cavity. The method then tests for convergence against a tolerance, and steps are iterated until convergence is achieved. The method then provides parameters of interest and the location of the cavity as output. Another aspect of this invention allows the calculation of cavity shape and cavity length for an input cavitation number. This is accomplished by an outer loop adjusting cavity length until the model converges to the input cavitation number. 

   
     BRIEF DESCRIPTION OF THE DRAWINGS 
     These and other features and advantages of the present invention will be better understood in view of the following description of the invention taken together with the drawings wherein: 
       FIG. 1  is a diagram of a partially cavitating axisymmetric body related to the method of the current invention; and 
       FIG. 2  is a flow chart of the method of the current invention; and 
       FIG. 3  is a flow chart of another method of the current invention. 
   

   DESCRIPTION OF THE PREFERRED EMBODIMENT 
     FIG. 1  shows a diagram of the physical problem of partial cavitation.  FIG. 1  shows a radial cross section of an axisymmetric body  10 . Axis r represents the radius from the axis of body  10 . Axis x represents the length along the body  10  measured from a cavitator disk  12 . Although a cavitator disk is shown, the model can calculate cavities for cavitator cones as well as cavitator disks. Flow, U ∞ , is in the direction of arrow  14 . A cavity  16  is shown extending from the edge of the cavitator along the length of body  10 . The length of the cavity, l c , is shown by dimension arrows. Likewise, the length of the body, l b , is also shown by dimension arrows. 
   Body  10  extends beyond a cavity closure  18 . Cavity  16  is closed to the body  10  with a modified Riabouchinsky cavity termination wall. Cavity closure  18  can be positioned in either body conical section  22  or body cylindrical section  24 . The plane of cavity closure  18  is referenced in the following disclosure as an endplate. 
   Body  10  has a flat front area  20  followed by a conical section  22  and a cylindrical section  24 . The diameter of flat front area  20  should be less than or equal to the diameter of the cavitator disk  12  base. 
   The mathematical formulations in of this algorithm are based on using the cavitator diameter to remove dimensionality for all lengths and using the free stream velocity, U ∞ , to remove dimensionality for all velocities. Alternate formulations using standard units can also be developed. 
   The flow field is governed by Laplace&#39;s equation,
 
∇ 2 Φ=0  (1)
 
where Φ is the total potential which is the sum of free stream potential, Φ ∞ , and disturbance potential, φ, giving:
 
Φ=φ ∞ +φ  (2)
 
   The free stream potential is the product of the velocity and the distance, x. Because the equation has been non-dimensionalized, the velocity is 1, and the free stream potential, φ ∞ , is x. The disturbance potential, φ, also obeys Laplace&#39;s equation, giving:
 
∇ 2 φ=0  (3)
 
   The disturbance potential satisfies Green&#39;s third identity, yielding a Fredholm integral equation of the second kind along the cavitator, cavity, endplate and body. Thus, at any point, x, on the body-cavity surface, the disturbance potential can be computed from: 
               2   ⁢     πϕ   ⁡     (   x   )         =       ∯   S     ⁢       [         ϕ   ⁡     (   x   )       ⁢     ∂     ∂   n       ⁢     G   ⁡     (     x   ;     x   ′       )         -       ∂     ∂   n       ⁢     ϕ   ⁡     (   x   )       ⁢     G   ⁡     (     x   ;     x   ′       )           ]     ⁢     ⅆ   S                 (   4   )             
 
where x′ are the points where the sources and dipoles are distributed under the boundary element model;
 
   S is the body-cavity surface; and 
   G(x,x′) is the Green function. 
   The Green function is further identified as: 
               G   ⁡     (     x   ,     x   ′       )       =     1          x   -     x   ′                      (   5   )             
 
   The dynamic condition on the cavity boundary is derived from Bernoulli&#39;s equation. Along the cavity surface, this can be written as: 
                 p   ∞     +       1   2     ⁢   ρ   ⁢           ⁢     U   ∞   2         =       p   c     +       1   2     ⁢   ρ   ⁢           ⁢     U   S   2                 (   6   )             
 
where P ∞  is the free stream ambient pressure;
 
   ρ is the free field fluid density; 
   p c  is the pressure inside the cavity; and 
   U s  is the flow velocity at the cavity surface.
 
The flow velocity at the cavity surface can be obtained from equation (6) giving: 
               U   S     =       1   +   σ               (   7   )             
 
where σ is the cavitation number which is defined as: 
             σ   =         p   ∞     -     p   c           1   2     ⁢   ρ   ⁢           ⁢     U   ∞   2                 (   8   )             
 
The kinetic boundary condition is that no flow crosses the body-cavity boundary, 
                 ∂   ϕ       ∂   n       =     -     n   x               (   9   )             
 
where n x  is the axisymmetric body free-stream velocity power. The no net flux condition, 
                 ∯   S     ⁢         ∂     ϕ   ⁡     (   x   )           ∂   n       ⁢     ⅆ   S         =   0           (   10   )             
 
is also required to make the problem a determinate system.
 
   Total drag is calculated by adding the drag coefficients. The pressure drag coefficient, C p , at {overscore (x)} is calculated as follows:
 
 C   p =1 −U ({overscore (x)}) 2   (11)
 
The pressure contribution to the drag coefficient may then be computed as: 
               C   dp     =       4   π     ⁢       ∯   S     ⁢       C   p     ⁢     n   x     ⁢     ⅆ   S                   (   12   )             
 
The viscous contribution to the drag coefficient along the wetted portions of the conical and cylindrical body areas is calculated using the International Towing Tank Conference equation given by Newman,  Marine Hydrodynamics , MIT Press, Cambridge, Mass. 1980, for the friction coefficient, C f , at {overscore (x)} is as follows: 
               C   f     =     0.075       (       log   10     ⁡     (       R   ⁡     (     x   _     )       -   2     )       )     2               (   13   )             
 
where R({overscore (x)}) is the local Reynolds number. The total viscous drag coefficient, C dv , is: 
               C   dv     =       4   π     ⁢       ∯   S     ⁢       C   f     ⁢     s   x     ⁢     ⅆ   S                   (   14   )             
 
The base drag coefficient, C db , which is the component of pressure drag associated with the base of the body is: 
                 C   db     =       0.029   ⁢       (     2   ⁢     b   base       )     3           C   dv           ,           (   15   )             
 
where b base  is the body radius at the base. The total drag coefficient is then given by
 
 C   d   =C   dp   +C   dv   +C   db   (16)
 
The panels are distributed along the cavitator, cavity, endplate, and cylindrical body section aft of the cavity, according to the partial floor method, known in the art. The partial floor method optimizes the number of panels in accordance with requirements for getting good convergence. Non-uniform panel spacing is used in many locations, in order to reduce the number of panels without reducing the accuracy of the solution.
 
   During iteration, the end plate height is determined by integrating the cavity surface back from its detachment point on the cavitator, and the number and distribution of panels along the endplate changes according to the changes in the endplate height. Smaller panels are required at highly non-linear flow locations, such as the region near the cavitator. Panel distribution in the wetted body area after cavity closure  18  changes to keep the aspect ratio of the neighboring panels between 0.5 and 2.0, in order to ensure good accuracy of the results. 
   In following the method of the current invention, first an initial cavity is defined. An arbitrary initial cavity can be chosen as a cone extending from the cavitator edge to an assumed endplate height of 0.2 or 0.3 is sufficient for most cases. In this discussion, the endplate height is measured as the radial offset from the body surface to the last point of the cavity. By applying equation (4) on all panels along the cavity body surface, S, a system of equations is obtained. This system is solved for the disturbance potentials, φ, along the wetted portions of the boundary and on the Riabouchinsky endplate; the normal derivative of the disturbance potential along the cavity boundary; and the cavitation number. 
   The kinetic boundary condition given in equation (9) is applied along cavitator, endplate, and aft body to update the cavity shape. In order to update the cavity, the program calculates how much each panel has to be rotated to satisfy the no flow condition. The program starts with the first panel at the cavitator and shifts the aft most point of the panel in the radial direction which satisfies the calculated rotation. The panel is rotated with the aft most point. The foremost point of the next panel is then shifted to the same radius as the previous aft most point. This process is continued until the panel adjacent to the endplate is undated. The endplate height is adjusted to the aft most point of the aft cavity panel. The iteration continues until the kinetic boundary condition converges to within a tolerance, giving the cavity shape. 
   From the converged disturbance potential along S, the disturbance velocity components can be calculated: 
               u   x     =           ∂   ϕ       ∂   x       ⁢           ⁢   and   ⁢           ⁢     u   r       =         ∂   ϕ       ∂   r       .               (   17   )             
 
   Referring now to  FIG. 2 , there is shown a flowchart of the current invention. In the input step  30 , geometric and other system parameter data including the estimated cavitation number, the estimated cavity length and the convergence criteria is read. The routine then distributes boundary element panels along the cavitator, cavity, endplate, body extension in the conical section, body extension in the horizontal section, and the aft body. The panels are distributed in order to reduce the number of panels and get an accurate result. In the initialize step  32 , the algorithm calculates the unit dipole and unit source functions and initializes matrices for the influence functions with known boundary values wherever applicable. The formulate equations step  34  formulates matrices for each panel using the disturbance potential equation (4) and no net flux condition given in equation (9). The solve equations step  36  solves the matrices created in the formulate equations step  34  in order to obtain the unknown disturbance potential along wetted body sources, normal distributions of disturbance potentials along cavity surfaces, and the cavitation number. The compute forces step  38  computes velocity components such as those in equation (17) and drag coefficients: including pressure drag, equation (12); viscous drag, equation (14); and base drag, equation (15) from the solved equations. In the update cavity step  40 , the cavity is updated from the computed forces using the kinetic boundary condition of equation (9). Convergence on cavity shape is checked in the converges decision step  42 . If the cavity is not converged, the initialize step  32  is executed to calculate influence functions for the updated cavity and next iteration thus begins. Once the cavity has converged, the compute parameters step  44  computes various output parameters of the converged solution which include pressure drag, viscous drag, base drag, total drag, cavitation number, cavity length, maximum cavity radius, length of cavity to maximum radius location. The output results step  46  then provides the location of the cavity written as coordinates and the cavity&#39;s disturbance potential, disturbance potential gradient, and pressure coefficient. 
   The basic algorithm enumerated above provides cavity shape and cavitation number based on an input cavity length. In order to obtain cavity shape and cavity length for an input cavitation number, the embodiment of  FIG. 3  adds an additional series of iterations. The user inputs a cavitation number and an assumed cavity length. This embodiment follows the previous embodiment in converging on a new cavitation number, σ, for the assumed cavity length. In step  48 , if the new cavitation number is within a tolerance of the given cavitation number, parameters are computed, step  44 , and the results are provided, step  46 . Otherwise the embodiment proceeds to step  50  wherein the algorithm determines the relationship between the new cavitation number, σ, and the given cavitation number. In step  52 , cavity length is increased by a predetermined amount if the calculated cavitation number is lower than the initial cavitation number, and in step  54  the cavity length is decreased by a predetermined amount if the calculated cavitation number is greater than the initial cavitation number. The routine loops back to the initialize step  32  and recalculates the cavitation number for the new cavity length. Operation continues until the calculated cavitation number falls within a tolerance of the initial cavitation number, the cavity length has converged, as tested in step  48 . 
   Using this invention; partial cavitation for high-speed underwater bodies can be analyzed. As disclosed, the invention can analyze axisymmetrical bodies using two cavitator shapes, a disk and a cone; however, the invention can easily be modified to analyze other axisymmetric cavitator shapes. As disclosed the inventive method can converge on cavity length or cavitation number. Total drag is calculated by adding the pressure drag, viscous drag and base drag. The invention can also be utilized for studying the effects of body aft radius, body cone angle and body cone angle starting at the cavity closure if the closure is on conical section  22 . This method provides new information concerning the physics of cavitation which can be used in the design of cavitating vehicles. 
   In light of the above, it is therefore understood that within the scope of the appended claims, the invention may be practiced otherwise than as specifically described.