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

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

OCR for page 93

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
SESSION 3
WAVY/FREE-SURFACE FLOW: PANEL METHODS 3

OCR for page 93

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

OCR for page 93

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
Nonlinear Ship Wave Calculations Using the RAPID Method
H.C.Raven (MARIN, The Netherlands)
ABSTRACT
The RAPID (RAised-Panel Iterative Dawson) approach for solving the fully nonlinear wave-resistance problem is described. This method, developed at MARIN, uses an iterative procedure based on a Rankine panel method similar to that of Dawson. No convergence problems are usually met in the speed and fullness range of practical ships, and even for cases showing extensive wave breaking in reality. The validity of the method regarding the occurrence of wave breaking and steep wave slopes is discussed. The modeling of the flow off a transom stern and its implementation in RAPID and in the corresponding linearized problem is studied. The nonlinear solution is found to be realistic and very accurate numerically. A number of validation studies shows that the RAPID results are very good in most respects, and often mean a surprisingly large improvement compared to linearized solutions. The principal remaining problems are resolution of short wave components, certain geometric complications, and the absence of a 3D model for wave breaking.
1. INTRODUCTION
Wave resistance is one of the most important resistance components of a ship, generally contributing 20 to 80 % of the total resistance of a ship sailing in still water. Experience has shown that this resistance component is quite sensitive to modifications in the design, and large reductions of the wave resistance can often be obtained without any important sacrifice in deadweight or cargo capacity. The effect of a bulbous bow is a well-known example. A capability to predict and minimize wave resistance in the design stage of a ship is, therefore, very desirable. Theoretical studies of the wave resistance problem have been carried out already for a century. Since Michell's paper of 1898, many studies of the linearized free-surface potential flow problem have been made, and several alternative solution methods and linearizations have been proposed; some of them very useful for a certain range of applications. Particularly in the early seventies much work has been done on this problem. In 1977, Dawson [1] proposed his well-known method, which in a number of comparative studies turned out to give realistic predictions, and to be more practical than most other proposals.
Much work has been done since on validation and improvement of Dawson 's method. There is a large body of literature on this method, and the mathematical background of the “engineering approach” has been more and more clarified. A code for Dawson's method is available at many institutes all over the world. Still, the actual use in commercial ship design work is quite limited. While at MARIN the code DAWSON [2] is being applied to 70 to 100 hull forms per year in commercial work and important successes have been obtained, most other versions of the method seem to be rarely applied.
As a matter of fact, the routine use at MARIN has made us aware of several shortcomings that make application more difficult. In particular, the resistance predicted by the code is not very reliable; for slender ships the estimate is usually fair, but for full hull forms the wave resistance is underestimated or even negative (a phenomenon explained in [3]). It is, therefore, not safe to optimize a hull form on the basis of the predicted wave resistance alone, as this sometimes even incorrectly indicates the order of merit of design variations. But at the same time the wave pattern and profile, the hull pressure distribution and flow direction have turned out to be realistic enough to be used for improving the design. Hull

OCR for page 93

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
form optimization using DAWSON is thus an art that requires insight in the assumptions made in the theory, knowledge of fluid dynamics and wave making, and much experience. Only few have collected enough experience to fruitfully use Dawson's method.
Although, therefore, the shortcomings of Dawson's method do not prohibit its very effective use, a next step has to be made. Not only a greater accuracy of the predictions is desired, but also we should like to include several aspects of the flow that cannot be taken into account in DAWSON. Specifically:
the neglect of nonlinear terms and transfer terms in the free surface boundary condition introduces inaccuracies;
the hull form above the undisturbed waterline cannot be taken into account in a linearized method; this eliminates the effects of bow flare, flat sterns, surface-piercing bulbous bows and so on;
doubts exist on the modeling of the flow off an immersed transom stern.
All these inaccuracies and limitations can in principle be removed by solving the fully nonlinear problem, rather than applying the slow-ship linearization underlying Dawson's method. At several institutes the development of a method for solving the fully nonlinear free surface flow problem is now in progress. Both panel or boundary integral methods, and field equation methods are used. The latter usually deal with viscous flows and almost invariably unsteady problems (or, at least, transient solution methods), and will not be discussed here.
The majority of methods for free surface flows around ships assumes a potential flow. Using Green's identity the problem is conventionally recast in a boundary integral form. There are many boundary integral methods addressing unsteady problems, e.g. [4,5]. Although in most cases the main purpose of the developments is to solve seakeeping or other actually unsteady problems, they may be applied to steady problems as well. One particular method to be mentioned here is that of Cao, Beck and Schultz [6] with its unconventional but very successful use of point sources at a distance above the free surface.
The most direct way to address the wave resistance problem is the solution of the steady free surface potential flow problem, which has been given much less attention. Most methods in this class are based on an approach similar to the linearized method of Dawson. This means that a distribution of Rankine sources or other simple singularities is used both on the hull and on a part of the free surface around it. For solving the nonlinear problem an iterative procedure is needed. Among the first to study such a method was Ni [7]. His method was Dawson-like, with panels on the free surface that were moved in every iteration. Although convergence was problematic, he succeeded to get converged solutions for a range of cases, subject to several restrictions on the use of particular difference schemes and so on. His method was later extended and improved in Larssons group [8]. Later Kim and Lucas [9] came up with a very similar approach, but found it necessary to apply artificial numerical damping in order to obtain convergence. In [10] they proposed a modified double iteration scheme which was claimed to improve the convergence. Another method was proposed by Delhommeau [11], also using panels on the free surface. Convergence again turned out to be problematic.
At the same time, some other developments were based on the use of singularity distributions not on but above the free surface. This choice was inspired by obvious practical advantages for the problem considered. One of the first efforts was due to Xia [12], who describes some initial attempts to solve the nonlinear problem. However, a too large sensitivity to the distance from the panels to the water surface made him give up. Much more successful were Jensen et al [13,14] with their method using point sources above the wave surface, analytical second derivatives and a special way of satisfying the radiation condition. In several further publications the method has been extended to shallow water, channel flows, unsteady cases in the frequency domain etc. At the same time, convergence of the method still can provide severe problems as soon as e.g. different panel distributions are used.
In [15,16], the present author proposed the so-called RAPID approach (RAised Panel Iterative Dawson). The use of free surface source panels raised above the free surface is largely similar to Jensen's method; but most details are different. The iterative procedure, selected on the basis of previous numerical estimates of the different nonlinear contributions [3], has appeared to be unusually stable and robust: Contrary to almost all previous methods for the same problem, no convergence problems are generally met. It thus seems to be one of the most practical methods, and, together with Jensen's method, it seems to be the only code for the nonlinear problem that is already actually applied in ship design work.
This paper describes recent progress in the further development and validation of the RAPID code. In [16], the essence of the method was explained and the results of several numerical experi-

OCR for page 93

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
ments were discussed. This will briefly be summarized in the next section, paying some more attention to the convergence properties and to the importance of nonlinear effects. Section 3 is devoted to a discussion of the flow off an immersed transom stern. The various flow regimes are studied, and the implementation of the transom flow model in both the linearized (DAWSON) and the nonlinear (RAPID) code are described and compared. Section 4 then discusses a number of detailed validation studies.
2. RAPID
2.1 The mathematical model
We consider the flow around a ship hull in a uniform flow from ahead. The origin of our coordinate system is located amidships at the level of the undisturbed free surface. The x-coordinate is positive astern, and y is positive in upward direction. The coordinate system does not move with the dynamic trim and sinkage of the ship. All quantities are nondimensionalized with the ship speed and length.
The potential flow assumption is made, and a velocity potential is introduced such that:
(1)
The field equation to be satisfied is the Laplace equation. The flow is undisturbed far ahead and at large lateral distances, while a wave pattern with a known rate of decay is present downstream.
On the wetted part of the hull the flow must be tangential to the hull, as represented by a Neumann boundary condition for the potential:
φn=0 (2)
On the free surface a kinematic boundary condition is to be satisfied, requiring that the flow is tangential to the wave surface:
φxηx+φyηy−φy=0 on y=η (3)
and a dynamic boundary condition that the pressure in the flow at the free surface must equal the atmospheric pressure (4):
on y=η (4)
2.2 Iterative solution
Contrary to what is done in linearized methods, these boundary conditions will be applied on the correct surfaces; i.e. the hull condition on the part of the hull that is actually wetted, and the free surface conditions on the actual wave surface. No simplifications are thus made. Solving this fully nonlinear free surface problem requires iteration. The iteration procedure is set up as follows.
While in a linearized method such as DAWSON the flow field and free surface are decomposed into a known base flow and a small perturbation, and linearization in terms of this perturbation is performed, we now make a similar assumption that the solution of a certain iteration differs only by a small amount from the flow field and free surface shape found in the previous iteration. Again we linearize in terms of this perturbation and impose the boundary condition at the previous free surface. Thus we obtain a linear fixed-domain problem to be solved in the current iteration. Having found the solution we update the base flow field and free surface, and start the next iteration.
Specifically, defining
η=H+η′ (5)
(6)
and substituting this into the dynamic and kinematic boundary conditions, we obtain linearized kinematic and dynamic conditions:
(7)
(8)
These two conditions are combined in exactly the same way as in Dawson 's method by eliminating the wave height perturbation, and we obtain the linearized free surface condition (FSC) to be imposed on the wave surface found in the previous iteration:
(9)
It is noted in passing that actually a consistent linearization would add certain terms representing the transfer of the boundary condition from the free surface to be calculated, towards the free surface found in the previous iteration. These terms are intentionally disregarded here, as they have been found to be rather irregular for the source distribution employed here and might lead to divergence of the iteration. To the converged solution these terms

OCR for page 93

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
make no contribution of course, since η' then vanishes and the FSC is imposed on the actual free surface.
Once the linear problem of the current iteration has been solved, the velocities found are substituted in the original dynamic condition (4), and a new free surface is found. This, together with the new velocity field, is the basis for the linearization of the next iteration.
This iterative process is started from any initial guess. As no assumptions have been made on the relation between the base flow field and the base free surface, this initial guess can be selected on the basis of convenience only. In most cases we start simply from a uniform flow and flat free surface; the first iteration thus becomes the Neumann-Kelvin problem. The iteration is repeated until the residual errors in both the dynamic and the kinematic boundary conditions are below a predefined tolerance in all free surface points (i.e. a maximum norm is used). In most cases, we impose tolerances of 0.1 % of the ship speed for the kinematic condition, and 0.2 % of the stagnation pressure for the dynamic condition. It is important to note that checking these residuals is the only safe procedure; visual inspection of the changes in a hull wave profile, as is often done, does not guarantee that the solution has converged !
2.3 Implementation
To implement the mathematical model outlined above, singularity distributions are specified for all boundaries where non-trivial boundary conditions are to be imposed, i.e. the wetted part of the hull and the surrounding part of the wave surface. On the hull the familiar distribution of quadrilateral constant-strength source panels is used, with collocation points in the panel centroids. For the free surface however, a less conventional choice has been made: quadrilateral constant-strength source panels are located at a certain distance above the wave surface, while the collocation points are on the wave surface itself. Within certain limits for the distance of the panel to the free surface this is permitted theoretically, and it has several practical advantages:
as there is freedom in the distance from a panel to the corresponding collocation point, the panels need not be repositioned in each iteration; only the collocation points are moved towards the new free surface position. This saves some work and avoids potentially destabilizing geometric manipulations;
the velocity field in the fluid domain is smoother than with constant-strength source panels on the boundary itself;
the numerical dispersion (i.e. the error in the wave length due to the discretization of the continuous velocity field and boundary condition) is far smaller than with the usual implementation;
the outer boundaries of the free surface domain are less reflective.
the implementation is simplified on several points.
There are, however, some disadvantages as well:
the conditioning of the system of equations is slightly worse; by a proper selection of the panel elevation this can, however, be minimized.
if the section shape has a large slope above the waterline, the collocation points right under the centers of the free surface panels adjacent to the hull are at a rather large distance from the waterline, with resulting inaccuracy. If the collocation points are kept at a constant distance from the actual waterline, the free-surface boundary condition is better resolved, but the conditioning is worse. This trade-off can present some problems, in particular for flat sterns, fter all, the use of raised panels has turned out to be extremely effective for this class of problems.
The limits of the distance from the panels to the free surface are, as shown in [16], roughly between 0.5 and 1.5 panel dimension. Here the panel length is more important than the transverse size. It will be obvious that the distance cannot be kept inside these limits if the panels are not adjusted to the free surface shape. E.g. the stagnation height at a ship speed U is U2/2g, and for a usual panel density of 20 panels per transverse wave length this is 1.6 times the panel length. Therefore, as described in [16], the free surface panel distribution must be adapted to the wave surface such as to stay at a more or less constant distance. But there is no need to do this very accurately, as the distance only has to be kept inside the limits mentioned before. Therefore, the calculation is made in separate jobs of 2 to 6 iterations, each time followed by an automatic adaptation of the free surface paneling to the last calculated wave shape, and of the hull paneling to the modified intersection with the free surface panels. If desired, the position of the hull is adapted to the new estimate of the trim and sinkage. Then the process is restarted from the results of the previous iteration. In this way, the efficiency and stability of the process is kept while a good conditioning and accuracy can be guaranteed.
Another detail of the implementation to be mentioned is the use of an upstream difference scheme for the derivatives of velocities in the free surface condition. This is exactly similar to the use in Dawson's method, and helps the satisfaction of the radiation condition at the expense of a small

OCR for page 93

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
numerical damping. The detrimental effect of this damping can be minimized by choosing a sufficient panel density. On the other hand, in order to increase the numerical accuracy at a given panel number, the use of a difference scheme might be dispensed with in the future.
The complete iterative procedure can be summarized as follows:
Define the initial hull and free surface panel distribution; choose the distance of the panels above the free surface.
Choose an initial free surface velocity distribution and wave elevation.
Impose the combined FSC, linearized with respect to the last free surface shape and velocity field, in the free surface collocation points; and the hull boundary condition in the hull points;
Solve the linear problem. Calculate the new free surface and velocity field;
Move the free surface collocation points to the new wave elevation, applying underrelaxation if desired;
If needed, adapt the free surface paneling to the new wave surface; adjust the hull paneling; adjust trim and sinkage of the hull.
Calculate the residual errors in all points, and decide to stop or to return to step 3.
2.4 Convergence and calculation time
As described in [16], the method displays
a good convergence with respect to increasing panel density;
only very little effect of the distance of the panels to the free surface;
a robust and easy convergence of the iterative procedure, without any tuning, smoothing or filtering.
Figure 1. Convergence of hull wave profile
Fig. 1 illustrates the convergence of the hull wave profile for a Series 60 Cb=0.60 hull at Fn=0.316, with a very dense paneling. 12 iterations were required to reduce the residuals to the desired level, but already after 4 iterations the hull wave profile is very close to the final solution.
An example of the convergence history is shown in Fig. 2, in which the logarithms of the maximum residual errors in the kinematic and dynamic boundary conditions are plotted against the iteration level.
Figure 2. Convergence history with and without free surface panel adaption
The convergence rate is constant after the initial iterations, and the residuals are reduced until a minimum level determined by the machine accuracy is reached. For the kinematic residual this level is higher than for the dynamic residual as can be understood from its more complicated expression.
The same case was run with the method without any free surface panel adaptation. From Fig. 2 the convergence appears to be slightly faster; the free surface panel adaptations in the original method at iterations 4, 8 and 13 slightly delay the error decrease as a result of e.g. the geometric changes. On the other hand, adaptation from time to time is unavoidable if one wants to use a fine paneling for steep waves and still have convergence. It is noted that methods using panels on the free surface itself are forced to apply such panel adaptations in every iteration. That a raised-panel method makes it unnecessary to do so may be one of the explanations for its robustness.
The code currently runs on a CRAY-YMP supercomputer. Calculation times vary widely according to the number of panels. To give some examples:
# hull panels
# FS panels
CP sec/iteration
340
540
2.7
680
1056
17
624
2312
100

OCR for page 93

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
These calculation times can probably be reduced substantially, but this has not been given priority up to now.
The number of iterations to reach the prescribed level of residual errors usually varies between 5 and 15, dependent on the difficulty of the case. On the basis of a visual inspection of the results, the calculations would probably be terminated after 5 iterations in most cases.
2.5 Nonlinearities and wave breaking
At the 1992 Symposium on Naval Hydrodynamics, discussers of [16] indicated some restrictions of this class of nonlinear solution methods. In the first place, methods that describe the free surface as a single-valued function of the two horizontal coordinates cannot cope with extreme wave slopes, which was stated to be an important drawback in practice. Thus such methods could impossibly converge in cases in which wave breaking occurs in reality, according to these discussers. And even, because of the above, this class of solution methods would be limited to so-called “weak singularities”. Similar statements can be found in [17].
As not all of these statements seem equally true, it may be useful to pay some more attention to the validity of the present method. Let us consider the physical phenomena that occur for steep waves. Most steady breakers are similar to spilling breakers: a region of turbulent flow with air entrainment, riding steadily on top of the wave pattern that otherwise remains relatively smooth. There are several experimental indications that this steady breaking occurs already at a smaller wave steepness than the Stokes limiting wave; but even the latter is described as a single-valued function of the horizontal coordinate. According to Duncan [18], wave breaking induced by a submerged hydrofoil occurs for wave slopes exceeding 16 to 24 degrees. That perhaps the spilling breaker itself has a very large local surface slope is a different issue, as the breaking will require another modeling anyway. But we can conclude that the potential flows giving rise to this kind of wave breaking and underlying the breaking zone can be easily described using a single-valued free surface; and that excessive wave slopes that cannot be described will rarely occur in reality because wave breaking eliminates them.
The next question is then: is it possible to compute waves so steep that breaking occurs? As mentioned above, some discussers expected that divergence would prohibit such calculations. In principle, I see three possibilities:
the iteration does not converge;
it converges, but a grid-independent solution cannot be obtained;
it converges, but to a non-breaking solution that is not physically realizable.
Figure 3 Calculated wave pattern for a tanker at Fn=0.25
Based on results obtained up to now the third option seems most likely, and we can at least reject the first possibility. As an example, Fig. 3 shows the wave pattern found for a tanker with a block coefficient of 0.83, at the unrealistic Froude number of 0.25. A converged solution was fairly easily obtained, although there is no doubt that extensive wave breaking will be present, e.g. in the wave trough behind the bow wave: The calculated local wave slope is almost 30 degrees, and the ratio of peak-to-trough wave height to wave length is 0.16. Both exceed known threshold values for inception of breaking. The panel length was 0.04 transverse wave length here, but grid independence has not been checked thoroughly for this case.
Similarly, in [10] a converged solution was obtained for a Series 60 CB=0.80 model at Fn=0.30, although on a rather coarse free surface grid.
Therefore wave breaking does not necessarily lead to divergence of the iteration process, and nonlinear methods like the one presented here can be a useful tool also in cases with wave breaking. But for completing the flow picture a separate criterion is required to determine whether a potential flow solution in practice will give rise to wave breaking, and some model of the breaking should eventually be incorporated. A proposal for such a model in 2D is [19].
The above example also illustrates that RAPID is stable enough to solve the nonlinear potential flow problem for all realistic speeds. At the same time it refutes the statement that these methods can only deal with “weak nonlinearities”. I am not aware of the existence of a clear definition of the latter, but I would not call the nonlinearities in this example

OCR for page 93

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
weak. The calculated peak-to-trough wave height is as large as the draft of the vessel and vertical velocities exceed 45 % of the ship speed.
Of course, other types of wave breaking occur. A sharp ship bow may act as a sort of ploughshare and induce a flow similar to a plunging breaker. This behaviour is definitely outside the reach of RAPID and all similar methods. Other methods, usually time-dependent, with a more general free surface description, have been applied to plunging breakers and do give a picture of the flow up to the moment of impact of the jet on the free surface. But then they stop ! Consequently also those methods have to be supplemented by some kind of local model for the plunging-type wave breaking and its effect on the flow field. In that respect they may have only a small advantage compared to the present approach.
An interesting related point of discussion is the question what happens if the discretization is further and further refined towards the bow. Again the same three pitfalls seem possible: divergence, unresolvable grid-dependence or an unphysical solution. The problem can be studied a little bit further by invoking the similarity between the steady flow at a sharp ship bow and a time-dependent 2D flow in cross planes. For a bow with a finite entrance angle the corresponding 2D case is that of an impulsively started wavemaker. From small-time expansions a logarithmic singularity in the wave elevation has been shown to exist in that case, representing a jet shooting up along the wavemaker. Numerical studies using time-dependent fully nonlinear boundary integral methods, e.g. [20,21] actually reproduced this behaviour. The finer the paneling, the higher the jet, and no grid-independence could be reached. But this is extremely localized, and the result at very small distances away from the wavemaker and everywhere else is totally unaffected by the grid dependence at the intersection. Analogously, and keeping in mind that in 3D steady flows singularities are likely to be weaker than in 2D time-dependent flows, I expect that for sharp bows in fact eventually no convergence will be obtained in RAPID if the discretization is refined towards the waterline; that locally no grid-independent solution or even no physically correct solution may thus be reached; but that the resulting error will be extremely localized and unimportant from the point of view of the global flow field. For almost all practical applications this theoretical drawback is then of no interest. Although there is some speculation in these statements, current experience with the method actually points in this direction. Further study is desired.
All this does not mean that the method proposed here is the final answer; but, subject to certain improvements in resolution and in numerical details, it may be the answer to most of the practical questions that can be dealt with using potential flow assumptions.
3. THE FLOW OFF A TRANSOM STERN
3.1 Physical Phenomena
As most current ship forms have a transom stern, the capability of a method to handle transom flows is of primary importance. In linearized methods the hull form above the undisturbed waterline has no effect, and transom sterns above the design waterline are usually not treated as transom sterns. But in a nonlinear method transom flow modeling is also required for all those cases in which a transom becomes immersed only at speed.
For ships having a transom stern, three flow regimes may be distinguished:
the “regular” type of flow, which is similar to the flow past a cruiser stern. This type of flow will occur if the transom is at a sufficient distance above the undisturbed free surface. Dependent on the hull form the last part of the actual waterline in potential flow may be either a streamline or an envelope of streamlines that leave the hull tangentially and continue on the free surface. At the smooth detachment of the free surface from the hull surface, both the hull boundary condition and the two free surface boundary conditions must be satisfied. Unless the detachment is in a stagnation point, this means that the flow leaves the hull tangentially.
the “transom flow”, occurring for immersed transoms at a sufficiently high speed. At the transom edge the hull streamlines leave the hull and continue on the free surface. The negative hydrodynamic pressure needed to make the pressure at the transom edge atmospheric is achieved by an upward curvature of the flow immediately aft of the transom edge, resulting in a more or less marked “rooster tail”.
The flow around a deeply immersed transom at low speed. At the sharp convex corner of the transom edge, viscous flow separation occurs and a dead-water region aft of the transom is formed.
The transitions between the various flow regimes are practically relevant, as a dead-water region is to be avoided, and a transom that is at a too high position as well. Let us further consider what phenomena determine these transitions. Suppose a regular flow type occurs, with the transom above the

OCR for page 93

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
water surface, and we increase the draft of the vessel. The detachment position is likely to move aft and at some moment reaches the transom edge. The streamline shape then further adjust so as to maintain a zero pressure at the edge. While in regular flow the curvature is bounded at detachment, this need not be true in the transom flow.
For a further increasing transom depth or for decreasing speed the hydrodynamic pressure coefficient required to balance the hydrostatic pressure at the transom edge increases in absolute value, and the upward curvature of the streamlines must increase. A very steep forward face of the wave aft of the transom results, eventually giving rise to a spilling breaker following the transom at a short distance. At still lower speed or greater draught this closes upon the transom and forms a dead-water zone.
From these preliminary considerations, we can list the requirements of the “ultimate” nonlinear potential flow model:
it should contain a model for at least the regular and the transom flow;
it should indicate the transition from a regular to a transom flow and reversely;
it should preferably indicate the transition from transom flow to dead-water flow.
If not all transitions can be incorporated in the code, one may try both flow regimes and check for anomalies. We can immediately note that point III will perhaps still be impossible to realize. If it is true that wave breaking determines this transition, a proper breaking criterion will be indispensable. Furthermore one should be aware of possibly large viscous effects on this transition. Moreover, point II requires the ability to compute the regular flow for extremely flat afterbodies; this is a very high demand for any iterative nonlinear method, as a small change of the wave elevation then results in a very large (and destabilizing) shift of the waterline.
Below, we shall consider the modeling of the transom flow in linearized and nonlinear models, and discuss the question of the transitions.
3.2 Mathematical Modeling
There seems to be some confusion about the type of modeling appropriate for transom flows. The boundary conditions to be imposed are fairly evident: on the hull the flow again must be tangential, and on the entire free surface it must be tangential and have a zero pressure. That there are two boundary conditions on the free surface corresponds with the additional degree of freedom there, the elevation of the free surface. At the transom edge however, the confluence of boundary conditions poses some problems. The two kinematic conditions can be reconciled if the flow leaves the hull tangentially. But as the edge location is fixed, the degree of freedom due to the wave elevation seems to be lost. How to impose the dynamic boundary condition?
This apparent overspecification of the problem reminds of the flow around a lifting body, in which circulation provides the additional degree of freedom. The analogy with a Kutta condition, which also imposes an extra condition on the flow at a sharp corner, has sometimes been invoked. Consequently, trailing vortex wakes or dipole sheets are introduced, e.g. in [22]. But the success of such an implementation could not be shown convincingly; in [22] the use of source panels on the transom instead of a trailing vortex wake led to almost identical results.
Another proposal for modeling transom flows is [23], which deals with slender hulls at very high speeds. Also here the necessity of introducing trailing longitudinal vorticity on the free surface is advocated. However, this vorticity appears not to result from the presence of a transom but from the use of a reverse image of the flow above the still water surface, as is appropriate in the infinite-Froude-number limit. The same trailing vorticity would in this model be required for hulls without transom. The slender-body formulation adopted causes this boundary vorticity to extend to downstream infinity with constant strength.
To study the meaning of trailing vorticity for transom flows, let us first consider the flow around a two-dimensional body with a transom stern. The flow domain is bounded by the hull and the free surface which originates from the transom edge. Therefore the domain is singly connected, and there cannot be any unknown circulation around the body.
In a three-dimensional case the situation is similar. Is there any reason why the three-dimensionality would suddenly introduce spanwise vorticity? If it would, trailing longitudinal vorticity would in fact also be present. However, this vorticity would start at the transom edge, i.e. on the free surface, so it would remain on the free surface; on the boundary of the fluid domain therefore. There is no cut in the domain, no additional degree of freedom, and the boundary vorticity has no physical effect on the flow field in the domain at all; it just affects the form of the Fredholm equation to be solved.
Therefore, in my opinion the use of a trailing vortex wake for transom flows, although perhaps

OCR for page 93

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
numerically helpful in some models, cannot be given any physical meaning. Rather than pointing out the analogy with the Kutta condition we might consider free streamline theory, as applied to e.g. the cavity or wake flow behind a flat plate at large incidence. Although the neglect of gravity may well affect some of the conclusions, free-streamline theory tells that the flow leaves the transom edge tangentially, but with infinite curvature. Even more illuminating is [24], in which for a flat-ship linearization in a 2D free-surface potential flow around a semi-infinite body the behaviour at the transom edge is derived. Again, the conclusion is that the slope of the streamline is continuous (tangential flow), but that the curvature tends to infinity at the edge. Specifically, the free surface shape aft of the transom is described as follows:
(10)
where ηtr is the height of the transom edge and the coefficient C2, which can be expressed in the hull geometry, is the degree of freedom required to satisfy the dynamic condition. Although the details of this result may be a consequence of the necessary simplifications and linearization, this study suggests that all transom conditions can be satisfied without any physical vorticity; that the additional degree of freedom is the streamline shape at the transom edge; and that a more or less singular behaviour at the edge may result.
3.3 Implementation
We shall now discuss the implementation of transom conditions. Because of the mentioned lack of general agreement on the proper modeling it seems useful to study the implementation in a linearized method as well. This has been used already for several years in the DAWSON code and provided good guidelines for the nonlinear implementation which we shall consider subsequently.
It is interesting to note that implementing transom conditions in a linearized method is basically inconsistent. Linearization assumes that the variation of the hull shape near the undisturbed waterline is small; i.e. a small curvature, and in slow-ship theory also a small slope of the sections and buttocks. Only then the exact position of the hull—free surface intersection has a negligible effect on the flow and linearization is permitted. However, at a transom edge the hull curvature is infinite, and we want to have a particular hull—free surface intersection. Nevertheless a quite reasonable representation of transom flows appears to be possible in linearized methods as well.
In the hull paneling the transom is left open. The free surface panel distribution is adapted to the transom by the addition of a few extra strips of panels aft of the transom. One might think of imposing transom conditions in points right at the transom edge; but with the usual constant-strength source panels on the hull the velocity field is singular here in the discretized formulation. Therefore transom conditions are to be imposed in the first collocation points aft of the transom.
In the linearized problem we can only impose one single FSC, a combination of the kinematic and dynamic FSC. The usual FSC is derived by expressing the term ηx in the kinematic condition in the wave elevations η in the point considered and a few upstream points using a difference scheme; and expressing these η's in the velocities through the dynamic condition. In this process any abritrary pressure added to the dynamic condition drops out; in other words, the combined form actually means that in the free surface collocation points the velocity must be parallel to the isobar surface passing through that point, but the pressure corresponding to that isobar surface is eliminated. The wave elevation is only retrieved afterwards from the dynamic FSC. Thus the ordinary form of the FSC makes it impossible to impose any condition on the wave elevation at or near the transom itself, and the free surface will most likely start somewhere at the transom face.
The solution of this difficulty is well known: the FSC must be modified so as to avoid elimination of the pressure level. This is easily achieved by imposing in the first collocation points the condition that the velocity must be directed parallel to a line from the transom edge to the wave elevation at the point, e.g.
(11)
The same transom condition can be derived by substituting into the usual FSC the velocity required for zero pressure at the transom edge.
Instead of the two-point formula for ηx, higher-order expressions may be used. These require that an assumption is made on the local behaviour at the edge. Two alternatives have been tried:
the behaviour according to Schmidt [24], (10)
a simple Taylor expansion in the transom edge point

OCR for page 93

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
(61)
The latter forms reduce the number of nested loops, with a resulting improvement in computational efficiency.
5.6
Field point on the panel
Consider first the case where the field point is at w=z0 in physical coordinates, corresponding to t=0 in parametric coordinates. In this case
(62)
and the logarithm can be decomposed in the form
(63)
If the last term is expanded in a Taylor series the formal result is
(64)
where
(65)
With the coefficients of (65) denoted by pn, successive powers of p can be expanded using the algorithm (25–26). Thus
(66)
where
(k=2,3 …N; j=k, …, N) (67)
and
(n=1,2, …, N–1) (68)
Substituting (66) for each term in the series of (64), and collecting the homogeneous coefficients, gives the desired result
(69)
where
q0=log z1 (70)
and
(j=1,2, …, N) (71)
The source integral can be evaluated using (69), with the result
(72)
The normal dipole integral requires a separate analysis. Evaluating the indicated derivative in (44), and using (66),
(73)
To evaluate the inverse of the series in (73), (27 –28) are used to find the coefficients , with , and . The resulting coefficients are
(74)
and
(n=1, 2, …, N) (75)
The normal dipole integral is then given as

OCR for page 93

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
(76)
Except for the singular case j+m+n=0 the integrals are elementary, with the result
(77)
where the integrals Iμ are defined by (44). The singular case occurs only for n=0, where the contribution from the term where m=0 and j= 0 is given by
(78)
Thus, when n=0, (77) applies in the modified form
(79)
Depending on the definition of ϕ0 (as a principal-value integral or as the limit where w→z0 on a specified side of the panel), the integral in (79) is either real, or differs from a real quantity by ±πi. The only contribution to the real part of (79) is in the latter case.
For the more general case where the field point is on the panel at w=z(t0), with t0 an arbitrary point in parametric space in the open interval (t1,t2), the results above can be extended using the binomial theorem. Hereafter a tilde is used to denote the shifted variable t̃̈=t−t0, or the coefficients in a series involving this variable. Thus the panel is defined in physical space by analogy with (17), in the form
(80)
where
(n=0,1,2, …, M) (81)
In place of (69) we write
(82)
where the coefficients q̄n are evaluated from the relations
(n=1, 2, …, M−1) (83)
(k=2, 3, …, N; j=k, …, N) (84)
(85)
(j=1, 2, …, N) (86)
The contribution to the source integral from the logarithmic singularity in (82) is
(87)
The contribution from the remaining terms in (82) is
(88)

OCR for page 93

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
where
(k=0,1,2, …, N) (89)
Thus the source integral is evaluated in the more general form
(90)
Following a similar generalization for the normal dipole integral, (73)–(76) are replaced by
(91)
where
(92)
(n=1,2, …, M−1) (93)
The singular case j=0 simplifies again, as a result of the identity
(94)
Thus
(95)
This equation can be simplified by extending the definition (46) of the coefficients ζn, with the final result
(96)
where
(k=1, 2, …, N) (97)
Here the coefficients ζn are defined for n= 0, −1, −2,… as in (46), but with the lower limit of the sum m=1−n.
As in the discussion following (79), the only contribution to the real part of the dipole potential from the singular integral which in (95–96) is from the residue.
6
System of Equations
We have discussed the B-spline defined singularity strengths and calculations of the influence functions in the above sections. In this section, we are going to discuss the implementations of the system of equations. The collocation method will be discussed first, and then the Galerkin method.
6.1
The Collocation Method
In this system of equations, the boundary condition is satisfied at collocation points, and the unknowns to be solved are the B-spline polygon vertices of the dipole moments. The number of

OCR for page 93

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
collocation points on each panel can be arbitrary subject to the condition that the system is not underdetermined.
Consider the discretized Green's formula (15). The influence of the dipole singularities on the collocation point i is
(98)
where j,k is the kth degree B-spline coefficient of the dipole moment on the jth panel. The moments of the influence functions can be evaluated as described in the last section, and we will define as the kth moment of the dipole influence function on collocation point i by panel j. Therefore, the dipole influence of panel j on collocation point i can be written as,
(99)
From section 3, j,k is a function of the B-spline polygon vertices, therefore, (99) can be further expanded as follows:
(100)
where ν,m is the mth B-spline vertex of the dipole moment, and the coefficients can be obtained by expanding the basis functions
(101)
Now the left-hand-side of the discretized Green's formula (10) can be rewritten as follows:
(102)
Hence, the left-hand-side of the system of equations has been established.
We can calculate the influence of the sources, the right-hand-side of the system of equations, in a similar way. Since the strengths of the sources are known, the calculations are relatively simple.
Starting with equation (13), the source strength of panel j, σj, can be expressed as . We next define the source influence function such that , where is the kth moment of the source influence function on the collocation point i by the panel j. Then the total source influence on the collocation point i will be,
(103)
Equation (103) gives the right-hand-side of the system of equations.
Therefore, the discretized Green's formula becomes,
(104)
Equation (104) is the complete formulation for non-lifting problems. The effect of the wake and the Kutta condition on the formulation of the lifting problem will be discussed later in this section.
Defining , and , we obtain the system of equations in the matrix form [A][Φν]= [RHS]. In this system, there are L+N unknowns, and L*Nc equations (Nc is number of collocation points per panel). Therefore, we need more than one collocation point per panel to avoid an underdetermined system. Increasing the number of collocation points per panel, Nc, increases the computational time but also improves the stability of the solution. Based on numerical experiments we have selected Nc =3. Using a cubic B-spline with three collocation points on each panel, the left-hand-side matrix [A] is,

OCR for page 93

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
(105)
and the unknowns are
(106)
This is an overdetermined system, which can be solved by a least squares approach. The solution of this least squares system initially was carried out by both the Householder method and the normal equation approach [2], and both methods gave the same solution. In the present method, the Householder method is used.
To solve a lifting problem, the influence of the wake must also be included. For a two-dimensional, steady lifting problem, the strength of the wake sheet dipole is equivalent to a point vortex at the trailing edge, and the strength of this vortex should satisfy the Kutta condition. As shown by Morino [5], the strength of this vortex Г is equal to the the potential jump at the trailing edge,
(107)
where + is the dipole moment at the trailing edge of upper side (suction side), and − is the dipole moment at the trailing edge of lower side (pressure side). The system of equations thus becomes,
(108)
where wi is the wake influence function, i.e. the contribution from the last term of equation (8) to the ith panel on the foil.
When applying B-splines to define the dipole moment, these two values are the dipole moments at the “end points” of the B-spline, and are therefore equal to the first and the last of the B-spline polygon vertices.
(109)
Therefore, the wake influence functions can be coupled into the left-hand-side of the least squares system. However, this will cause the system of equations to become more ill-conditioned, and will therefore increase the numerical error in the solution. In the present method, we first guess the initial solution Г0, and then impose a pressure Kutta condition to get the final solution iteratively. The pressure Kutta condition used here was first introduced by Kerwin, et al. [3] for threedimensional propeller problems, and it requires the equality of the pressures on the upper and lower sides of the trailing edge. This non-linear pressure Kutta condition has been described in [ 3], and it can be simplified for two-dimensional problems. It is found that by using this approach, the solution procedure is more stable and more accurate. In the present paper, we used the two-dimensional flat plate solution (2πα, where α is the angle of attack) to be our initial solution, however, the selection of the initial solution (including zero) will not affect the final solution since the satisfaction of the Kutta condition should provide a unique solution.
For non-lifting problems, the accuracy is improved by imposing a continuity condition at the two ends of the B-splines. That is,
(110)
Although solving a least squares system sometimes is less accurate and robust than solving a square linear system, it provides more freedom to set up the system of equations. The selection of the number of collocation points on each panel is one example, and the applications of the Kutta condition and the continuity condition described above further demonstrate the flexibility of the present method. However, the size of the matrix is increased by the number of collocation points. This issue will be discussed later.
6.2
The Galerkin Method
The Galerkin method used here is the Bubnow-Galerkin method, in which the test functions are selected to be the same as the interpolation functions.

OCR for page 93

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
Applying the Galerkin method to the integral equation expressed in equations (8) and (9), we will have,
(111)
where the dipole moment, (t), can be expressed in terms of the B-spline basis functions as equation (7).
The inner integrals in the equation (111) are exactly the same as the collocation method described in the last section, therefore, the discretized form of equation (111) can be derived similarly. Let us define functions , and as the discretized forms of the inner integrals of the left-hand-side and the right-hand-side of the equation (111), where li stands for the lth collocation point of panel i. Therefore, they represent the total influence on the collocation point l of the panel i of the dipole and source (plus the wake effect for the lifting problems) separately. From equation (108), we have,
(112)
Equation (111) becomes the following equation by replacing these two functions with the inner integrals:
(113)
By applying a quadrature rule to this equation, we then have,
(114)
where ωl is the weighting function of the numerical integration. The resultant system of equations is an L+N by L+N linear system for a degree N higher order panel method. If we expressed the matrix system as equation (105–106), the matrix system of the Galerkin method will then be [Q][A][Φν]=[Q][B], and the matrix [Q] is,
(115)
if the number of collocation points is selected as 3.
The Galerkin method utilized here uses the exact form for the inner integrals, and a numerical integration for the outer integrals. Chandler and Sloan [1] called this a “qualocation” method.
The numerical examples shown later show that the Galerkin method has almost the same accuracy and efficiency as the collocation method. However, the left-hand-side matrix size of the Galerkin method will be smaller than that of the collocation method. For a cubic higher order panel method, with L panels and Nc collocation points, the collocation methods will give an L*Nc by L+3 left-hand-side matrix, and the Galerkin method will give an L+3 by L+3 left-hand-side matrix. If we used a Householder least squares solution for the collocation method, then it needs 2(L+3)2(L*Nc−(L+3)/3) flops. On the other hand, the Galerkin method needs 2(L+3)3/3 flops to solve the matrix system by using Gaussian elimination. The collocation method needs about 3(Nc−1) times more computational time than the Galerkin method. This factor is not important for two-dimensional problems because only a small portion of the computational time is spent on solving the matrix (less than 10%).
In the next section, we will show results for both the non-lifting and lifting problems by applying the present method.
7
Computational Results
In this section, several numerical examples will be shown for the present method. We will first show results of the collocation method for all the cases, and then compare the solutions of the Galerkin method to the solutions of the collocation method.

OCR for page 93

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
Figure 4: The number of panels used on one quadrant vs. the maximum error of the solutions of different order panel methods for the case of a uniform inflow past a square.
Figure 5: The absolute error (maximum error of the solutions) of the case a uniform inflow past a square, vs. the computational time of different order panel methods. The computational time is measured on a SiliconGraphics Iris 4D25 machine.
The first case we tested is that of uniform flow past a square. Since the present method is a general higher order method, we will present the results of linear potentials, quadratic potentials, and cubic potentials, and compare these results to the analytical results and results of a low order panel method.
Due to the sharp corners of the square, a non-uniform B-spline curve is used to define the geometry. An end condition is applied for the continuity as described in (110), such that the first B-spline control point coincides with the last B-spline control point. The error of this test case is defined as the maximum difference between the computational results and the analytical solution. Figure 4 shows the results of linear potential, quadratic potential, and cubic potential from the present method, along with the results from a two-dimensional constant potential panel method. The accuracy of the higher order panel method is much better than that of the low order panel method. Figure 5 shows the errors vs. the computational time (run on a SiliconGraphics Iris 4D25 machine) by using the cubic potentials and constant potentials. For a 1% error, the computational time of the constant potential method is about 36 times longer than that of the cubic potential method. This means that the present method is not only accurate, but also efficient.
The second case we tested is a symmetrical Karman-Treftz foil with a thickness/chord ratio of 25% and a tail angle of 45º (as shown in figure 1). The angle of attack is selected as 5 degrees. The error in this case is defined as the difference between the computational circulation and the analytical solution. Figure 6 shows the errors by using the different order of panel methods. It shows the similar trend as the square case. One of the advantages of applying the B-splines to define the solutions is that the derivatives of solutions can be directly calculated from the B-spline basis functions and solved B-spline polygon vertices. That is, the velocities can be computed directly without using finite differencing schemes. Figure 7 shows the calculated pressure distribution of this foil by using cubic potentials compared to the analytical results. It shows that even 10 panels gives a very good result for the pressure distribution. Notice that although we only show the pressure coefficients on the panel boundaries, in fact the solutions are continuous because they are defined by B-splines. Figure 8 again shows the error vs. the computational time of this case, and for the 1% error, the computational time of

OCR for page 93

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
Figure 6: The number of panels used vs. the error of the solutionsof different order panel methods for the Karman-Treftz foil case.
Figure 7: The pressure distribution of calculated results and analytical solutions of the Karman-Treftz foil case.
Figure 8: The absolute error (error of the calculated circulation) vs. the computational time of different order panel methods for the hydrofoil case. The exact solution is 0.33. The computational time is measured on a SiliconGraphics Iris 4D25 machine.
the constant panel method is 41 times that of the cubic potential method. Notice that the constant potential method needs 480 panels to get a result within 1% error, and the cubic potential method only needs 10 panels.
To explore alternative solution techniques we have compared the results based on the collocation and Galerkin methods. Three collocation points are used on each panel, which were selected initially to provide uniform spacing in the parametric coordinate t, at the points 1/6, 1/2, 5/6 of the distance between each pair of vertices. The Galerkin scheme was implemented using the same integration nodes with equal weights. Subsequently, to refine the Galerkin integration, Gauss nodes were substituted in both methods. As shown in Figure 9, the improved accuracy of the Galerkin method is practically insignificant, while the use of Gauss nodes reduces the error in both methods by a factor of order two.
Figure 10 demonstrates the error em vs. the number of elements in the LHS matrix (hereafter refer to as the “matrix size”) of the constant potential method, and higher order (cubic) collocation method and higher order Galerkin method.

OCR for page 93

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
Figure 9: The error of computational results from the collocation method and the Galerkin method. The test case is the Karman-Treftz foil
Figure 10: The error vs. the matrix size of different order panel methods for the Karman-Treftz foil case.
The matrix size of the constant potential method is the square of the number of panels, L, the matrix size of the higher order collocation method is 3*L (for three collocation points each panel), multiplied by L+3, and the matrix size of the higher order Galerkin method is the square of L+3. From figure 10, one can see that 480 panels of the constant potential, 10 panels of the collocation and Galerkin higher order methods have similar errors (about 0.7% of the exact value). However, the number of elements in the matrix of these methods are 230400, 390 and 169 separately. That is, for the same error, the matrix size that the constant potential method needs is about 600 times larger than that of the higher order collocation method, and 1400 times larger than that of the higher order Galerkin method! Therefore, the presented method also has the advantage of using less computer resources.
8
Conclusion
We have presented a higher order panel method by using B-splines to define both the geometry and solutions. The test solutions show that for a given number of panels, the present method is more accurate than the low order method, and for a given accuracy, the computational time is far less than that of a low order panel method. The solutions are naturally continuous, and the derivatives of potentials can be directly calculated from the B-spline polygon vertices and B-spline basis functions. Also this method needs far less computer storage space or memory than the low order method, which will be significant advantage when solving three dimensional problems.
This method has been implemented by both the collocation method and the Galerkin method. The accuracies of both methods are equally good. The left-hand-side matrix of the Galerkin method is smaller than the collocation method. The resulting savings of the computational time, memory, and storage space are relative minor in solving two-dimensional problems; however, these savings are expected to be substantial in solving three-dimensional problems.
The extension of this method to solve three dimensional problems is under research. The problems of three-dimensional geometries are much more complicated, such as the connection of two surfaces, or the intersection of two surfaces. These complexities have to be well defined, and carefully incorporated into the solution

OCR for page 93

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
matrix system. The calculations of the influence functions are also much more difficult than the two-dimensional problems. It is hoped that this method will not only be able to provide more accurate solutions than the low order panel method, but also be more effective of solving problems such as wave body interactions, and lifting problems.
Acknowledgements
This work was performed as part of the Joint Industry Project “Wave effects on offshore structures”, and also under the MIT Sea Grant College Program with support provided by the David Taylor Model Basin. The authors would like to thank Mr. Seamus Tuohy, a graduate student in MIT Ocean Engineering Design Laboratory, for his help in answering our questions about B-splines.
References
[1] G.A.Chandler and I.H.Sloan. Spline qualocation methods for boundary integral equations. Numer. Math., 58, pp. 537–567, 1990.
[2] G.H.Golub and C.F.Van Loan. Matrix Computations. The Johns Hopkins University Press, Baltimore, Maryland, 1989.
[3] J.E.Kerwin, S.A.Kinnas, J-T Lee, and W-Z Shih. A surface panel method for the hydrodynamic analysis of ducted propellers . Trans. SNAME, 95, 1987.
[4] Sir Horace Lamb. Hydrodynamics. Cambridge University Press, sixth edition, 1932.
[5] Luigi Morino and Ching-Chiang Kuo. Subsonic potential aerodynamic for complex configurations: a general theory. AIAA Journal, vol 12, no 2, pp. 191–197, February 1974.
[6] M.B.Okan and S.M.Umpleby. The use of B-splines for the calculation of two-dimensional potential around arbitrary bodies. Int. Ship-building Progress, 32, June 1985.
[7] D.F.Rogers and J.A.Adams. Mathematical Elements for Computer Graphics. McGraw-Hill, Inc., New York, 1990.

OCR for page 93

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