By: Robert L. Daily

Except where reference is made to the work of others, the work described in this thesis is my own or was done in collaboration with my advisory committee. This thesis does not include proprietary or classified information.

This study attempts to define an optimal hull shape for a boat towing either a skier or a wakeboarder. Two methods for determining a free surface deformation (both hull and wake shape) given an applied pressure disturbance are compared. The methods derive from the same potential flow, but vary in how the pressure disturbance is defined; one uses a Fourier type approximation, the other a piecewise constant interpolation.

Using the results from the comparison of the two methods, an approach for simulating the wake shape given different hull shapes (via the pressure distribution) is developed. Using this approach the hull shape is optimized based on two wake parameters: height and slope. Examples of this optimization are presented for different towing scenarios. The resulting hulls and wakes are then examined for realism.

v

**ACKNOWLEDGEMENTS**

The author would like to thank his advisor Dr. Peter Jones for his guidance through the highs and lows of this research. He would also like to thank the other committee members, Dr. Jay Khodadadi and Dr. A. J. Meir, for their instruction both in classes and when problems arose during research, not to mention their patience through this process.

The author additionally wishes to express his gratitude to his entire family for their support throughout his entire school career and for instilling a deep love for the activities motivating this research. Also, he wishes to thank Ms. Susan Wooden for her support and encouragement.

vi

Style manual or journal used: ASME Journal of Dynamics Systems, Measurement, and Control.

Computer software used: Microsoft Word 2003.

vii

TABLE OF CONTENTS

List of Figures………………………………………………………………………………………………………x

List of Tables…………………………………………………………………………………………………….xiii

Nomenclature…………………………………………………………………………………………………….xiv

1. Introduction…………………………………………………………………………………………………..1

2. Mathematical Derivations……………………………………………………………………………….7

2.1. Basic Assumptions and Relationships……………………………………………………….7

2.1.1. Conservation of Mass………………………………………………………………………9

2.1.2. Kinematic Boundary Condition……………………………………………………….10

2.1.3. Dynamic Boundary Condition…………………………………………………………11

2.1.4. Far Field Boundary Condition…………………………………………………………11

2.2. Velocity Potential………………………………………………………………………………….12

2.3. Sinusoidal Pressure Patches……………………………………………………………………14

2.3.1. Surface Integral……………………………………………………………………………..15

2.3.2. Wave Number Integral……………………………………………………………………21

2.3.3. Wave Angle Integral………………………………………………………………………31

2.3.4. Evaluation Grid……………………………………………………………………………..41

2.4. Constant Pressure Patches………………………………………………………………………45

2.4.1. Surface Integral……………………………………………………………………………..47

2.4.2. Wave Number Integral……………………………………………………………………48

viii

2.4.3. Wave Angle Integral………………………………………………………………………48

2.4.4. Evaluation Grid……………………………………………………………………………..49

2.5. Numerical Issues…………………………………………………………………………………..51

2.5.1. Sinusoidal Pressure Method…………………………………………………………….52

2.5.2. Constant Pressure Method………………………………………………………………59

3. Method Validation……………………………………………………………………………………….63

3.1. Sinusoidal Pressure Basis Function…………………………………………………………63

3.2. Simplified Realistic Hull Pressure…………………………………………………………..67

3.3. Sinusoidal Pressure Centerlines………………………………………………………………69

4. Method Issues……………………………………………………………………………………………..73

4.1. Constant Pressure Method Symmetry………………………………………………………73

4.2. Sinusoidal Pressure Method Symmetry……………………………………………………75

4.3. Sinusoidal Pressure Method Oscillation…………………………………………………..77

4.4. Issues With Inversion…………………………………………………………………………….82

5. Hull Optimization………………………………………………………………………………………..88

5.1. Event Definition……………………………………………………………………………………88

5.2. Typical Hull Shape……………………………………………………………………………….92

5.3. Test Cases……………………………………………………………………………………………94

5.4. Jump 1 Test Results………………………………………………………………………………97

5.4.1. Resulting Pressure and Hull Shape…………………………………………………..98

5.4.2. Resulting Wake……………………………………………………………………………100

5.5. Jump 2 Test Results…………………………………………………………………………….102

5.5.1. Resulting Pressure and Hull Shape…………………………………………………103

ix

5.5.2. Resulting Wake……………………………………………………………………………105

5.6. Slalom 1 Test Results…………………………………………………………………………..106

5.6.1. Resulting Pressure and Hull Shape…………………………………………………107

5.6.2. Resulting Wake……………………………………………………………………………108

5.7. Slalom 2 Test Results…………………………………………………………………………..109

5.7.1. Resulting Pressure and Hull Shape…………………………………………………110

5.7.2. Resulting Wake……………………………………………………………………………111

5.8. Wakeboard 1 Test Results……………………………………………………………………112

5.8.1. Resulting Pressure and Hull Shape…………………………………………………113

5.8.2. Resulting Non-Optimal Pressure and Hull Shape……………………………..114

5.8.3. Resulting Non-Optimal Wake………………………………………………………..115

5.9. Wakeboard 2 Test Results……………………………………………………………………116

5.9.1. Resulting Pressure and Hull Shape…………………………………………………117

5.9.2. Resulting Wake……………………………………………………………………………118

5.10. Results Summary…………………………………………………………………………….119

6. Conclusions and Suggestions for Future Research………………………………………….121

6.1. Conclusions………………………………………………………………………………………..121

6.2. Suggestions for Future Research…………………………………………………………..123

References………………………………………………………………………………………………………..126

x

LIST OF FIGURES

Fig. 1 Nautique Ski 196…………………………………………………………………………………………2

Fig. 2 Monterey 268 SS SuperSport………………………………………………………………………..3

Fig. 3 Super Air Nautique 210………………………………………………………………………………..4

Fig. 4 Coordinate Frame and Uniform Flow…………………………………………………………….8

Fig. 5 Sinusoidal Pressure Patch Definition……………………………………………………………15

Fig. 6 Complex Plane and Line Integrals………………………………………………………………..23

Fig. 7 Sinusoidal Pressure Hull Discretization………………………………………………………..42

Fig. 8 Constant Pressure Hull Discretization…………………………………………………………..46

Fig. 9 Sinusoidal Pressure Combined Integrands…………………………………………………….55

Fig. 10 Sinusoidal Pressure II Integrand with Slow Oscillations……………………………….56

Fig. 11 Sinusoidal Pressure Final Integrands…………………………………………………………..58

Fig. 12 Constant Pressure Original Far Field Integrand……………………………………………60

Fig. 13 Constant Pressure Far Field Integrand with Slow Oscillations……………………….61

Fig. 14 Constant Pressure Far Field Final Integrand………………………………………………..62

Fig. 15 Sinusoidal Pressure Basis Function Test Pressure………………………………………..64

Fig. 16 Sinusoidal Pressure Basis Function Test Wake…………………………………………….64

Fig. 17 Sinusoidal Pressure Basis Function Test Hull………………………………………………65

Fig. 18 Simplified Realistic Hull Test Pressure……………………………………………………….67

Fig. 19 Simplified Realistic Hull Test Wake…………………………………………………………..68

xi

Fig. 20 Simplified Realistic Hull Test Hull…………………………………………………………….69

Fig. 21 Sinusoidal Pressure Basis Function Centerlines……………………………………………71

Fig. 22 Single Constant Pressure Patch Wake…………………………………………………………74

Fig. 23 Corresponding Constant Pressure Patch Wakes……………………………………………74

Fig. 24 Total Constant Pressure Patch Basis Function……………………………………………..75

Fig. 25 Sinusoidal Pressure Patch Wakes……………………………………………………………….76

Fig. 26 Loss of Data as More Fourier Terms are Used……………………………………………..78

Fig. 27 Oversampled Pressure Basis Function………………………………………………………..79

Fig. 28 Oversampled Wake Basis Function…………………………………………………………….80

Fig. 29 Oversampled Wake Basis Function Surface………………………………………………..81

Fig. 30 Original Test Surface………………………………………………………………………………..83

Fig. 31 Pressure Resulting From Original Test Surface……………………………………………84

Fig. 32 Pressure Resulting From Rotated Test Surface…………………………………………….85

Fig. 33 Increasing Pressure Oscillation………………………………………………………………….87

Fig. 34 Ski Jumping Example and Course Layout…………………………………………………..89

Fig. 35 Slaloming Example and Course Layout………………………………………………………90

Fig. 36 Wakeboarding Examples…………………………………………………………………………..91

Fig. 37 Wake Parameters……………………………………………………………………………………..92

Fig. 38 Basic Planing Hull Shape………………………………………………………………………….93

Fig. 39 Examples of Planing Boats………………………………………………………………………..95

Fig. 40 Jump 1 Minimization………………………………………………………………………………..98

Fig. 41 Jump 1 Pressure and Hull Shape………………………………………………………………100

Fig. 42 Jump 1 Wake…………………………………………………………………………………………102

xii

Fig. 43 Jump 2 Minimization………………………………………………………………………………103

Fig. 44 Jump 2 Pressure and Hull Shape………………………………………………………………104

Fig. 45 MasterCraft X-Star…………………………………………………………………………………105

Fig. 46 Jump 2 Wake…………………………………………………………………………………………106

Fig. 47 Slalom 1 Minimization……………………………………………………………………………107

Fig. 48 Slalom 1 Pressure and Hull Shape…………………………………………………………….108

Fig. 49 Slalom 1 Wake……………………………………………………………………………………….109

Fig. 50 Slalom 2 Minimization……………………………………………………………………………110

Fig. 51 Slalom 2 Pressure and Hull Shape…………………………………………………………….111

Fig. 52 Slalom 2 Wake……………………………………………………………………………………….112

Fig. 53 Wakeboard 1 Maximization…………………………………………………………………….113

Fig. 54 Wakeboard 1 Pressure and Hull Shape………………………………………………………114

Fig. 55 Wakeboard 1 Non-Optimal Pressure and Hull Shape………………………………….115

Fig. 56 Wakeboard 1 Non-Optimal Wake…………………………………………………………….116

Fig. 57 Wakeboard 2 Maximization…………………………………………………………………….117

Fig. 58 Wakeboard 2 Resulting Pressure and Hull Shape……………………………………….118

Fig. 59 Wakeboard 2 Wake………………………………………………………………………………..119

xiii

LIST OF TABLES

Table 1 Test Parameters……………………………………………………………………………………….94

Table 2 Typical Tow Boat Parameters……………………………………………………………………96

NOMENCLATURE

A amplitude coefficient of a pressure patch (N/m2)

B half–beam (m)

E leading edge contribution tensor (-)

*E final wave height basis tensor (-)

I total number of longitudinal solution points in the wetted area (-)

J total number of transverse solution points in the wetted area (-)

L wetted length of the craft (m)

M total number of longitudinal pressure patches (-)

N total number of frequency multiples / transverse pressure patches (-)

P pressure (N/m2)

U vehicle speed (m/s)

a pressure patch length (m)

b pressure patch width (m)

g acceleration due to gravity (m/s2)

i longitudinal index of solution point (-)

j transverse index of solution point (-)

k solution point index (-)

k0 fundamental wave number (1/m)

m longitudinal index of pressure patch (-)

xiv

n frequency multiple of the pressure patch / transverse index of pressure patch (-)

p pressure strip index (-)

t time (s)

u x coordinate of velocity (m/s)

v y coordinate of velocity (m/s)

w z coordinate of velocity (m/s)

x forward coordinate (m)

y port coordinate (m)

z up coordinate (m)

ˆi forward direction (-)

ˆj port direction (-)

ˆk up direction (-) v total fluid velocity (m/s)

Φ total flow potential (m2/s)

ζ height of the fluid surface (m)

η coordinate of the pressure (-) ˆj

θ wave angle (rad)

ξ coordinate of the pressure (-) ˆi

ρ fluid density (kg/m3)

φ craft disturbance potential (m2/s)

xv

1

1. INTRODUCTION

Wake shape is an important factor to consider when designing recreational or professional water skiing towboats. With the gain in popularity of wakeboards and other water sports, the “ideal” wake has become much more complicated. Most professional ski events are based around objects in the water: slalom buoys, jump ramp, etc [AWSA, IWSF]. For these events the competitor wants as little interference from the wake as possible, i.e. a flat wake. Wakeboard events, on the other hand, rely on the wake for the competition. Most events are based around using the wake as a launching platform for varying jumps [WWA, WWC]. For these events the competitors want a wake shaped to allow higher jumping. The shape of the wake is affected by several boat characteristics and operating conditions such as trim, speed, and mean draft. The largest design variable is the hull shape.

Most current towboats are designed based on experience and past designs. The potential for radical improvements to boat design exists in utilizing computer simulations to test designs before they are manufactured. There are two philosophies to designing towboats; the first is pure design (either wakeboard or ski). These boats are designed for the sole purpose of towing a specific piece of equipment. Other considerations such as passenger room and comfort are secondary to towing. The boat is designed to create the best wake for the specific equipment. A boat with this design is usually reserved for tournament tow boats; they are too limited for general public enjoyment. An example of

this type of design is the Nautique Ski 196. Note the closed bow and limited seating. This boat is not designed for passengers. Also notice the engine compartment in the center of the boat; this configuration is known as an inboard engine. Directly in front of the engine compartment is the ski pylon where the tow rope is attached. The forces from the rope act near the center of the boat and so pull it from its course less. It is not obvious, but the hull is also designed specifically with the wake in mind.

Fig. 1 Nautique Ski 196

The second philosophy is mixed design. These boats are for more general recreational use. The wake produced can be used for either skiing or wakeboarding, but is not ideal for either. An example of this type of design is the Monterey 268 SS SuperSport. This boat has an open bow and much more seating. The engine compartment is at the rear of the boat in what is known as an inboard/outboard configuration. There is no ski pylon interior to the boat to attach the tow rope; instead on the transom there is a hook. Because the rope attach point is at the rear of the boat, when the skier swings wide it can cause the rear of the boat to also move in that direction. The hull of this boat is designed more for passenger comfort than wake performance.

2

Fig. 2 Monterey 268 SS SuperSport

Wakes shapes for this type of boat are very hard to optimize. The main differences in operating conditions between skiing and wakeboarding are rope length (rider placement in the wake) and boat speed. The goal of the design is to make the wake act differently based on these two conditions. Recently some companies have started putting water tanks at the stern of boats such as the Super Air Nautique 210; these tanks increase both draft and trim when full allowing the boat to create a bigger wake for wakeboarding. Another feature of most pure wakeboarding tow boats is the tower which takes the place of the ski pylon. The tow rope attaches much higher allowing for better jumps, etc.

3

Fig. 3 Super Air Nautique 210

A design method using computer simulations would significantly reduce the cost and increase the performance of the hull design. The traditional method to simulate fluid flow is computational fluid dynamics (CFD). This method involves numerically solving the Navier-Stokes equations [Currie]:

()()()2Ptρρμλμ∂+−∇−∇+∇=+∇∇∂uuuuugg ρf. (1.1)

with μ and λ viscosity parameters of the fluid. In general this is not an easy problem to even simulate. One simplification that can be made is to assume inviscid and incompressible flow resulting in Euler’s equations:

()Ptρρ ρ

∂++∇=∇∂uuug f . (1.2)

This simplification removes some terms from the differential equations, but they are still nonlinear and without a closed-form solution. Several methods for numerically approximating the solution to this equation are available. These methods work very well for fluids in enclosed spaces. They provide in-depth knowledge of both the pressure and the velocity of the fluid throughout the entire volume in question. Determining the wake

4

5

shape does not require that knowledge, however; it only requires the free surface height. The height can be determined from the fluid pressure and velocity information but several problems with the CFD methods exist. One problem is the boundary conditions: to numerically approximate the solution the conditions at the boundary of the volume in question must be known. In enclosed spaces this is generally a condition based on the velocity of the fluid at a solid wall or a velocity or pressure profile at an inlet or outlet. In open water this is not as easy to determine. Related to the boundary conditions is a bigger problem, the free surface itself. One boundary of the fluid is the free surface. CFD methods require the boundary location to be known. However, for the wake shape problem the air-fluid boundary location is the unknown. To solve this problem using CFD an initial free surface height must be guessed and then the solution iterated by varying the shape and position of the free surface until all known boundary conditions are met. This would take an inordinate amount of time, particularly because the free surface consists of many different points each of which must be varied independently. In conclusion, a different method needs to be employed to determine the wake shape a boat produces.

Researchers have studied the flow around planing bodies for many years. The initial research centered on determining equations for the flow by using different simplifications for the hull shape [Sottorf, Maruo 1951, Maruo 1967]. This research was used to calculate the loading on the hull and thus the horsepower requirements for the boat [Savitsky]. The research evolved to using numerical technique approximations instead of simplifying the hull [Wang, Shen, Doctors]. Two numerical approaches evolved; the first uses a series of vortices to represent the flow [Lai]. The more common approach is to use

6

pressure elements to represent the hull on the surface [Tong, Cheng]. This idea has proven successful for determining lift and drag for hulls [Wellicome]; it is only recently being used to also determine the free surface height behind the boat [Scullen].

The idea to determine the wake shape is to take a pressure distribution on a free surface with a known wake shape. If this wake shape is linear in the pressure distribution then many of the pressures can be combined over the wetted area of the boat. The pressures are varied so the deformation of the surface in the wetted area matches the hull shape. The free surface height behind the wetted area from the resulting pressures is then the wake shape. Ideally, the method could be inverted so that a given hull shape (or perhaps even desired wake shape) produces the pressure. This research compares two pressure functions used as the basis for the pressure distribution. It then extends the method into finding the wake created by using many pressure patches. The method presented allows any pressure to be modeled and therefore any hull should be possible. Using the developed method an optimal hull shape for both skiing and wakeboarding given the constraints each activity is under is examined.

2. MATHEMATICAL DERIVATIONS

As the introduction mentioned, traditional computational fluid mechanics methods are not appropriate for computing the wake a planing craft creates. Instead, a technique that creates basis functions linear in the pressure amplitude will be utilized. These basis functions are summed together with different pressures to create the hull shape of the craft. The pressure that created this hull is then used to create the wake shape.

2.1. BASIC ASSUMPTIONS AND RELATIONSHIPS

To create the surface basis functions, the following assumptions are made about the fluid and craft:

1. The craft is traveling at a constant speed (U). This assumption is equivalent to a stationary craft in a flow with constant velocity.

2. The fluid is inviscid.

3. The fluid is incompressible.

4. The fluid has infinite depth.

5. The disturbance created by the craft is small compared to the constant velocity flow.

The flow created under these conditions can be represented as a potential flow. The potential of the total flow (Φ) is given by:

UxφΦ=−. (2.1) 7

where φ is the potential of the disturbance created by the craft. x () is the direction the craft is traveling, y () is to the port side of the craft, and z () is up. ˆiˆjˆk

Fig. 4 Coordinate Frame and Uniform Flow

y

x

U

From the definition of a potential flow the total fluid velocity is determined by

ˆˆˆuvw=++=∇Φvijk. (2.2)

The following notation will be used to denote derivatives throughout the paper.

2,xxyxxy∂Θ∂Θ=Θ=

Θ

∂∂∂ (2.3)

where Θ is a generic scalar variable. Along with these basic assumptions, several preliminary properties are required to derive the velocity potential.

8

2.1.1. Conservation of Mass

The first property to examine for the flow is conservation of mass [Currie],

0VDdVDtρ=∫. (2.4)

ρ is the fluid density, t is time, and V is the volume of interest. Expanding this equation using Reynolds’ Transport Theorem,

()tVVDdVdVDtΘ=Θ+∇Θ∫∫v, (2.5)

the conservation of mass equation (2.4) then becomes

()()()0txyzVuvwdVρρρρ⎡ ⎤ +++=⎣∫ ⎦ . (2.6)

This equation must be true for an arbitrary volume, so the integrand must always be zero,

()()()0txyzuvwρρρρ+++=. (2.7)

Expanding the derivatives in this equation gives

()0txyzxyzuvwuuuρρρρρ++++++=. (2.8)

Because the fluid is assumed incompressible (0αρ=) the conservation of mass equation is further simplified to

0xyzuvw++=. (2.9)

In terms of the perturbation potential, the conservation of mass equation is

0xxyyzzφφφ++=. (2.10) 9

2.1.2. Kinematic Boundary Condition

The next property to examine is the kinematic boundary condition. It states that the normal velocity of the surface (DzDt) is equal to the change in height of the surface with respect to time (DDtζ) where ζ is the height of the fluid surface. This idea is written as

()0DzDtζ−=. (2.11)

Expanding this using the definition of a material derivative,

txyDuvwDt y

Θ=Θ+Θ+Θ+Θ, (2.12)

the kinematic boundary condition (2.11) becomes

()()()()(txyDzzuzvzwzDt

)z

ζζζζ−=−+−+−+−ζ . (2.13)

The assumption that the craft is traveling at a constant speed implies the system is in steady state,

()0xyzuvw1 ζζζ=−−+−. (2.14)

Written in terms of the perturbation potential,

0xxxyyzzU z φζζφζφζ=−+−+− φ

z

. (2.15)

Because of the assumption that the perturbation potential (and therefore the change in surface height) is small, it is safe to linearize this equation. Additionally, it is assumed that this equation applies on the undisturbed free surface (z = 0) instead of the true surface (z = ζ). Therefore, the kinematic boundary condition is

0xUζφ=+. (2.16) 10

2.1.3. Dynamic Boundary Condition

The next boundary condition to examine involves the conservation of momentum (Bernoulli’s equation) applied on the surface of the fluid. The Bernoulli equation states

12PGFρ+∇Φ⋅∇Φ−= (2.17)

where F is the Bernoulli constant, G is the gravity term (Ggz =−), and P is the pressure. On the free surface the Bernoulli equation is,

(2222122xyzxPFU U ) g φφφφζρ=+++−++. (2.18)

To determine the Bernoulli constant, the Bernoulli equation is evaluated far upstream where the surface is undisturbed (ζ = 0) the pressure is atmospheric (P = 0) and the velocity is constant (ˆU=−v i). The resulting Bernoulli constant is

212FU=. (2.19)

Substituting this constant into the Bernoulli equation (2.18) and linearizing results in the dynamic boundary condition:

0xPUgφζρ=−+. (2.20)

This equation will be used to determine the wave height, so the last term is left in. In all other aspects it is assumed to apply on the undisturbed free surface, however.

2.1.4. Far Field Boundary Condition

The final boundary condition to consider is the far field boundary condition. This condition derives from the assumptions about the fluid far away from the disturbance created by the craft. In particular, the condition holds in two locations. The first is far

11

below the craft (z=−∞); the second is far ahead of the craft (x=∞). The condition states that the velocity perturbation will not affect the vertical velocity of the fluid ( 0 zφ=).

2.2. VELOCITY POTENTIAL

The conditions presented thus far apply to any velocity potential that meets the assumptions presented in section 2.1. Wehausen and Laitone derive a velocity potential for regular waves that meets these assumptions and describes the fluid behavior created by applying a pressure to the free surface of the fluid. This potential is

()()()()()22202023sec0220011,sec2secsincoscossinseccosseccossecsinkzSkzkePUkkkxkydkdkekxkydddπππθπφξηθπρπθξθηθθθξθηθθθξη∞−−⎧⎪=⎨−⎪⎩−−−⎡⎤⎡⎤⎣⎦⎣⎦⎫⎪⎡⎤−−⎡⎤⎬⎣⎦⎣⎦⎪⎭∫∫∫∫∫gg (2.21)

where S is the wetted area of the applied pressure, ξ and η are the coordinates of the pressure, θ is the wave angle, k is the wave number, and k0 is the fundamental wave number defined as

02gkU=. (2.22)

To determine the surface height, the dynamic boundary condition (2.20) is used. This equation requires the partial derivative of the velocity potential with respect to the longitudinal direction on the free surface. This partial derivative is

12

()()()()()22200202240220011,2seccoscoscossinsecsinseccossecsin.xzSkPUkkkxkydkdkkxkydddππππφξηπρπθξθηθθ θ

ξθηθθθ∞=−−⎧⎪=⎨−⎪⎩−−+⎡⎤⎡⎤⎣⎦⎣⎦⎫⎪⎡⎤−−⎡⎤⎬⎣⎦⎣⎦⎪⎭∫∫∫∫∫gg ξ η

(2.23)

Solving the dynamic boundary condition (2.20) for surface height results in

()()()()()()()2220202240220011,,2seccoscoscossinsec,sinseccossecsin.SkxyPgkkkxkydkdkPxykxkydddgππππζξηπρπθξθηθθθξθηθθθξηρ∞−−⎧⎪=⎨−⎪⎩−−+⎡⎤⎡⎤⎣⎦⎣⎦⎫⎪⎡⎤−−−⎡⎤⎬⎣⎦⎣⎦⎪⎭∫∫∫∫∫gg (2.24)

An equation for wave height has now been developed; all that remains is to define the pressure acting on the surface and evaluate the integrals.

The double integral over the surface does not have an analytical solution for all pressures. The solution could be solved numerically at this point, but that would involve solving four nested integrals numerically for each point in a solution grid. This is very expensive computationally. To reduce the computation time pressures are chosen that result in an analytical solution for the double integral over the surface. Two forms of pressure functions will be presented. Both are linear in the pressure amplitude coefficients. The first method creates a pressure that varies sinusoidally in the transverse direction (y, η) and is piecewise constant in the longitudinal direction (x, ξ). The second method is piecewise constant in both the longitudinal and transverse directions.

13

2.3. SINUSOIDAL PRESSURE PATCHES

The sinusoidally varying pressure method utilizes pressure patches that vary as sine waves in the transverse direction to approximate the true pressure [Cheng]. Multiple pressure patches over the same surface with different sine frequencies are added together. This is similar to a Fourier approximation of a function. The total pressure is piecewise constant in the longitudinal direction. In general, the form of the pressure is

()()1(,)sin2NnnnPABBBπξηξηη=⎡⎤=+−⎢⎥⎣⎦Σ ≤ ≤ B. (2.25)

The frequency of the sine wave is modified by changing n; An(ξ) is the amplitude coefficient of the pressure patch, B is the half-beam, and L is the wetted length of the craft. With this pressure, the surface equation (2.24) becomes

()()()()()()()011011,,,2;1sin,21,,,;2BNnnnBLNnnBNnnnBLASxyddgByBnAxyBxygByBASxyddBygξξηξηπρπζρξξηξηπρ=−−==−−⎧⎪⎪−≤≤⎪⎪⎡⎤−+=⎨⎢⎥⎣⎦⎪⎪≤−⎪≤⎪⎩Σ∫∫ΣΣ∫∫ (2.26)

with

()()()()()()222020224002201,,,sincoscos2seccossinsecsinseccossecsin.nnkSxyBkxBkkkydkdkkxkydπππππξηηξπθ

θ

ηθθθξθηθθθ∞−−⎧⎪⎡⎤=+−⎡⎤⎨⎣⎦⎢⎥−⎣⎦⎪⎩−+−⎡⎤⎡⎣⎦⎣⎫⎪⎡⎤−⎬⎣⎦⎪⎭∫∫∫gg

⎤⎦

(2.27)

14

Fig. 5 Sinusoidal Pressure Patch Definition

x = -L

x = 0

x = -a

m = 1

m = 2

m = M

y = B

y = 0

y,η

x,ξ

a

The first step in determining the wave height is to evaluate the double integral over the wetted area.

2.3.1. Surface Integral

The pressure is piecewise constant in the longitudinal direction, i.e. An(ξ) is piecewise constant in ξ:

15

a

()(),1,2,,1,0;0;0;2;1;2;0;nnnnmnMnMAaAaaAAamamALaLaALLLξξξξξξξξ−<⎧⎪−<<⎪⎪−<<−⎪⎪⎪=−<<−−⎨⎪⎪−+<<−+⎪⎪−<<−+⎪⎪<−⎩MM. (2.28)

The length of a pressure patch is given by a, and m defines the particular longitudinal patch. Therefore, the surface height becomes

()()()()()()()()()0,1,21212,,11,1,,,,,,2,,,,,,sin2,,,.BaNnnnnnBaaamLanmnnMnamLaNLannnMnLxyASxydASxydgASxydASxynAxyBBASxyddgζξηξπρ

,

d

ξ η ξ

ξηξξηξπξηξηρ−=−−−−−−+−−−+−+=−⎡=+⎢⎣++++⎡⎤+⎢⎥⎤⎣⎦+−⎥⎦Σ∫∫∫∫∫Σ∫LL (2.29)

To simplify notation, it is assumed the transverse location is in the half breadth of the craft (ByB−≤≤). If this is not the case the last term of the equation drops out. The condition will be added back later. An intermediate variable is used to evaluate the longitudinal integration:

()*1amξξ=+−. (2.30)

Therefore, the generic longitudinal integral appears as

()()()(10*,,,,,1,amnnamaSxydSxyamd ) * ξηξξηξ−−−−=−−∫∫. (2.31)

16

The integrand of this integral (2.27) is

()()()()()()()()()}()()()22*2020*224*00220***,1,,1,sin2scos1coscossinsecsin1seccossecsin,,1,,,,.nnnmnkSxyamBBkkkxamkydkdkkxamkydSxyamSxyπππππξηηπθξθηθθ

ec

θξθηθθθξηξη∞−−⎧⎪⎡⎤−−=+⎨⎢⎥−⎣⎦⎪⎩⎡⎤−+−⎣⎦−⎡⎤⎣⎦⎡⎤+−+−⎣⎦⎡⎤−⎣⎦−−=∫∫∫ggg (2.32)

Substituting this integrand back into the height equation (2.29) results in

()()()()()()()()()00******,1,1,2,2100******,,,1,10***1,,,1,,,,,,,2,,,,,,sin2,,,1,2BNnnnnnBaanmnmnMnMaaNnnnMnManmxyASxydASxydgASxydASxydnAxyBBASxyddgxyAgζξηξπρξηξξηξπξηξηρζπρ=−−−−−−−=−⎡=+⎢⎣++++⎡⎤+⎢⎥⎤⎣⎦+−⎥⎦=Σ∫∫∫∫∫Σ∫LL()()

ξ η ξ

()0***1,11sin2,,,.NBnNMnnmnmBanAxyBBSxyddgπξηξηρ===−−⎡⎤+⎢⎥⎣⎦−ΣΣΣ∫∫(2.33)

The surface integrand (2.32) is split into two terms, the single and the double integral:

()()()******,,,1,,2,,,,,,,,,nmnmnmSxySxySxyξηξη=+ ξ η . (2.34)

17

with

()()()()()()()()()()22**,,12020*2**24,,202*2001,,,sin2seccos1coscossin,,,sinsec2sin1seccossecsin.nmnmnkSxyBBkkkxamkydkdnSxykBBkxamkydπππππξηηπθξθηθπξηηθ

θ

ξθηθθ∞−−⎡⎤=+⎢⎥−⎣⎦⎡⎤−+−−⎡⎤⎣⎦⎣⎦⎡⎤=+⎢⎥⎣⎦⎡⎤⎡⎤−+−−⎣⎦⎣⎦∫∫∫

θ

(2.35)

This implies the surface integral term of the wave height equation (2.33) is also split into two parts:

()()()()(,,1,,2111,,,sin2NMNnmnmnnmnngxygxygxyAxyBBπρζρζρζ===⎡⎤⎡⎤=+−⎣⎦⎢⎥⎣⎦ΣΣΣ + ) (2.36)

with

()()()()()()()()()()0202,,,122020**220,4,,22*020,sin22scos1coscossin,sinsec22sin1seccossecsBnmnmBaBnmnmBaAnkxyBBkkkxamkydkdddAknxyBBkxamkyπππππζη

ec πθξθηθθξηπζηπξθηθ∞−−−−−−⎡⎤=+⎢⎥−⎣⎦⎡⎤−+−⎣⎦−⎡⎤⎣⎦⎡⎤=+⎢⎥⎣⎦⎡⎤−+−⎣⎦−∫∫∫∫∫∫∫gggg*in.dddθθξη⎡⎤⎣⎦

θ

(2.37)

These two equations are integrated with respect to the pressure coordinates (ξ and η) resulting in

18

()()()()()()()()()()2,,,1220201,11,22,12,22,03,,221,11,22,12,2,sec4sec,,,,,sec4,,,,nmnmnmnmAkxykkIkIkIkIkdkdAkxyJkJkJkJkdππππζθπθθθθθζθπ

θ

θθθθ∞−−=−⎡⎤+−−⎣⎦=⎡⎤+−−⎣⎦∫∫∫gg θ

(2.38)

with

()()()()()()()()()()()()()(),2200,20121sinsin,1sin1,221,21cosseccossec,1secsin21cossincossin.nllljjnllljjkkIknklBjkkJknkBxamyBxamyBωωθπθωθωθθπθθωθθωθθ+−+−±±−−=−+==−−=−+=+−+⎡⎤⎣⎦=++mm (2.39)

It is important to note that when the pressure strip is rotated by the wave angle θ, 1ω± and 2ω± correspond to the leading and trailing edges of the strip respectively. The wave height component equations (2.38) are written in terms of the leading and trailing edge contributions as

()()()()()(),,1,,1,,,1,,,2,,2,,,2,,,,,nmnmLnmTnmnmLnmT

,

,

xyxyx y

xyxyζζζζζζ=−=− x y

(2.40)

19

with

()()()()()()()()()()2,,,1,1,11,2220202,,,1,2,12,2220202,03,,2,1,11,22,0,,2,,sec,,4sec,sec,,4sec,sec,,4,nmnmLnmnmTnmnmLnmnmTAkxyIkIkkAk

k dkd

xyIkIkkAkxyJkJkdAkxyππππππ

k dkd

ζθθπθ

θ θ

ζθθπθζθθθθπζ∞−∞−−⎡⎤=+⎣⎦−⎡⎤=+⎣⎦−⎡⎤=+⎣⎦=∫∫∫∫∫

θ θ

()()232,12,22sec,,.4JkJkdππθθθθπ−⎡⎤+⎣⎦∫ (2.41)

The trailing edge is equivalent to the leading edge with the longitudinal coordinate (x) shifted by the pressure patch length (a).

()()21, , xyxaωω±±=− y, (2.42)

which implies

()()()(),,1,,,1,,,2,,,2,,,,,nmTnmLnmTnmL .

xyxa y

xyxaζζζζ=−=− y

(2.43)

Therefore, each of the wave height components (2.40) only depends on the leading edge contribution:

()()()()()(),,1,,1,,,1,,,2,,2,,,2,,,,,nmnmLnmLnmnmLnmL

,

, .

xyxyxa y

xyxyxaζζζζζζ=−−=−− y

(2.44)

Overall, the total surface height is the sum of the heights made by each pressure patch and the static deflection from the pressure being applied, i.e.

()()()(,111,,sin2NMNnmnnmnngxygxyAxyBBπρζρζ===

) ⎡⎤=− + ⎢⎥⎣⎦ΣΣΣ. (2.45)

20

The contribution of each pressure patch is written as a leading edge component and a trailing edge component. The trailing edge component is identical to, but offset from, the leading edge component:

()()(),,,,,,,nmnmLnmLgxygxygxayρζρζρζ=− − , . (2.46)

The leading edge component is split into two terms, a single and a double integral:

()()(),,,,1,,,2,,,nmLnmLnmL , . xyxyζζζ=+ x y (2.47)

The wetted area integral is now evaluated. Two integrals remain: the wave number and wave angle integrals. Only the double integral term contains the wave number integral. To make the wave angle integration similar for both terms, the wave number integral is evaluated next.

2.3.2. Wave Number Integral

To evaluate the wave number integral in the double integral term of the wave height (2.41) a complex integration method must be employed. The basic form of the integral in this method is

1,2001,2secjjkEIdkkkθ∞=−∫ j = . (2.48)

With I1,j from equation (2.39) this integral becomes

()()()()()2001sinsin1,2sec1sin2nmmjjkkkEnkkkBωωπθθ+−∞⎡⎤−−⎣⎦==⎡⎤−−+⎢⎥⎣⎦∫ dk j . (2.49)

21

To evaluate this integral a complex variable and equation are defined:

()()()1121iQjjQkiuQeFQQksQkω±±=+=−−−− (2.50)

where s is simply a variable substitution and k1 and k2 are values at which the function becomes singular. To match the integrand,

2102sec(1)2sinsin.jkknkBsθπθθ==−−=− (2.51)

On the real axis this function evaluates to

()()()1121ikjjkeFkkkskkω±±=−−−−. (2.52)

It can also be written in terms of sine and cosine functions,

()()()()()112cossin1jjkikkFkkkskk

1 ωω±±±+=−−−−, (2.53)

instead of an exponential. This equation is beginning to look like the integrand presented above. In fact the integral in equation (2.49) is

()()()()()11001sgnImsgnImnjjEFkdkFkdkωω∞∞++−−

j

⎧⎫⎧=−−

⎫

⎨⎬⎨⎩⎭⎩⎭∫ ∫ ⎬

j

. (2.54)

This integral is divided into two parts:

()1njjEee+−=−− (2.55) 22

with

23

k dk

()()10sgnImjjeFω∞±±±⎧⎫=⎨⎬⎩⎭∫. (2.56)

In the complex plane, the line from zero to infinity is divided into five parts. Part one lies on the real axis and goes from the origin to immediately before the first singularity. Part two circles the first singularity with an infinitely small radius. Part three is also on the real axis between the two singularities. Part four circles the second singularity similar to part two. Part five is on the real axis again and goes from the second singularity to infinity. Adding two more parts extends this line to a closed contour. Part six is a quarter circle from the real axis to the imaginary axis with an infinite radius, and part seven lies on the imaginary axis from infinity back to the origin.

u

Fig. 6 Complex Plane and Line Integrals

Q – plane

k1

l6

l5

l3

l2

l1

R→∞

k2

l4

r→0

r→0

k

l7

The real integral from zero to infinity is equivalent to an integral over segments one, three, and five:

24

j Q dQ

()()1351sgnImjllleFω±±+++⎧⎫⎪⎪=⎨⎬⎪⎪⎩⎭∫. (2.57)

Using the Cauchy Theorem, the line integral over a closed contour is zero,

()0cjlFQ±=∫ (2.58)

and the integral over the real segments written in terms of the other segments is

()()1352467jlllllllFQFQ±+++++=−∫∫ j

±

0

. (2.59)

Therefore, equation (2.57) becomes

. (2.60) ()()()()()24671sgnImjjjjjlllleFQdQFQdQFQdQFQdQω±±±±±±⎧⎫⎪⎪=−+++⎨⎬⎪⎪⎩⎭∫∫∫∫

To evaluate this total integral, each of the component integrals is studied separately.

The first integral to examine is over the l2 segment. The following parameters define the line to integrate over:

. (2.61) 21:,0iilQkredQiredrββββπ=+=≤≤→

The integral therefore becomes

()()()()1120101112lim1iikreijijiirlkreeFQdQiredkrekskrekβωββββπβ±+±→+=−−+−+−∫∫. (2.62)

Evaluating the limit, gives

()()()11201121ikjjlikeFQdQdskkωπβ±±=−−−∫∫, (2.63)

or

()()()1121121ikjjlikeFQdQskkωπ±±−=−∫. (2.64)

In terms of sine and cosine functions the integral is

()()()()()211111121sincosjjlkFQdQkikskkπω±±−

ω ± ⎡⎤=−+⎣⎦−∫. (2.65)

Taking the imaginary component of this integral leads to

()()()(2111121ImcosjjlkFQdQkskkπ ) ω±±⎧⎫−⎪⎪=⎨⎬−⎪⎪⎩⎭∫. (2.66)

The next integral to consider is over the l4 segment. The following parameters define this segment:

42:,0iilQkredQiredrββββπ=+=≤≤→0. (2.67)

The integral over this segment is:

()()()()2140202122lim1iikreijijiirlkreeFQdQiredkrekskrekβωββββπβ±+±→+=−−+−+−∫∫. (2.68)

This integral is split into two conditions, one where the singularity is positive and one where it is negative. If the singularity is positive (k2 > 0), the limit evaluates to the l2 component:

()()()21402121ikjjlikeFQdQdskkωπβ±±−−=−∫∫. (2.69) 25

This integral is solved as

()()()2142121ikjjlikeFQdQskkωπ±±−=−∫. (2.70)

Written in terms of sine and cosine functions the integral is

()()()()()422121121sincosjjlkFQdQkikskkπω±±−

ω ± ⎡⎤=−+⎣⎦−∫. (2.71)

If the singularity is negative (k2 < 0), the limit evaluates to

()40jlFQdQ±=∫. (2.72)

Therefore, the total imaginary component of this integral is

()()()()422121221cos;0Im0;jjlkkkFQdQskkkπω±±⎧−⎧⎫>⎪⎪⎪=−⎨⎬⎨⎪⎪⎪⎩⎭ 0 <⎩∫. (2.73)

To simplify evaluation of this function, the inequality condition is converted into terms of wave angle over which this function will be integrated. From equation (2.51) the singularity being positive (k2 > 0) implies

(1)02sin(1)0.sinjjnBπθθ−−>−< (2.74)

j is given, so this function solved explicitly in terms of the wave angle is

1;022;0.2jjπθπθ=<<=−<< (2.75) 26

The singularity being negative (k2 < 0) implies

(1)02sin(1)0.sinjjnBπθθ−−<−> (2.76)

This function is also solved explicitly in terms of the wave angle as

1;022;0.2jjπθπθ=−<<=<< (2.77)

The next integral to consider is the l6 integral. The following parameters define the line to integrate over:

6:,02iilQRedQiRedRββπββ==≤≤ → ∞. (2.78)

The integral therefore becomes

()()()162102lim1iiReijijiiRlReeFQdQiRedReksRekβωπβββββ±±→∞=−−−−∫∫. (2.79)

Written in terms of sine and cosine functions this integral is

()()()16sincos22120lim1RiijjRiilieeFQdQdkkseeRRωββπββββ±−+±→∞=−−⎛⎞⎛⎞−−⎜⎟⎜⎟⎝⎠⎝⎠∫∫. (2.80)

Much of the integrand is simplified by taking the limit, but one term remains:

()()()162sincos01limRijjRliFQdQedsπωβββ±−+±→∞=−−∫∫. (2.81)

27

This term depends on the value of 1ω±. If 10ω±= the integral is solved as

()()612jjliFQdQsπ±−−=∫. (2.82)

On the other hand, if the integral evaluates to 10ω±≠

()60jlFQdQ±=∫. (2.83)

Therefore, taking the total imaginary component of this integral leads to

()()6111;Im20;jjlFQdQsπωω±±±⎧−−⎧⎫⎪⎪⎪ 0

0

==⎨⎬⎨⎪⎪⎪⎩⎭≠⎩∫. (2.84)

Once again, it is simpler to write the condition in terms of the wave angle. From equation (2.39), this condition becomes

()()()()()11cossin01cos1tan.xamyBxamyBxamyB

sin

ωθθθθθ±=+−+⎡⎤⎣⎦=+−+⎡⎤⎣⎦−+−=mmm (2.85)

Solving this condition explicitly for the wave angle produces

()11tanxamyBθ−−+−⎛⎞=⎜⎝⎠m ⎟ . (2.86)

The last integral to consider is over the l7 segment. The following parameters define this segment:

7:,0lQidQidγγγ==≤ < ∞. (2.87)

28

This integral is

()()()170121jjleFQdQdiksikγωγγγγ±−±∞=−−−∫∫. (2.88)

Through some manipulation, this is written as

()()()()()()()()()()()117121112122222121201122222212120111.jjljkkeeFQdQdskkkkkkkkekkeidskkkkkkγωγωγωγωγγγγγγγγ±±±±−−∞±−−∞⎡⎤⎢⎥=−−−+−−⎢⎥++⎣⎦⎡⎤⎢⎥−−−−−⎢⎥++⎣⎦∫∫∫ (2.89)

The variable substitution ikγλ= is made producing the following integrals:

()()()()11111222001222001sgnsgn.1iiiikikiiieeddgkkkeedkdkfkkγωλωγωλωγλγλωγλγλγλ±±±±−−∞∞±−−∞∞±==++==++∫∫∫∫ iω

(2.90)

These solutions involve the even and odd auxiliary functions to the sine and cosine functions [Abramowitz]:

()()2020;01;01ztzttegzdtztefzdtzt∞−∞−

.

=≥+=≥+∫∫ (2.91)

Using these integrals, the imaginary portion of the total l7 integral is

()()()()()()()()71111122221121Imsgnsgn.jjlkFQdQkfkskkkkfkkkωω±±±⎧⎫⎡−−⎪⎪=⎨⎬⎢−⎪⎪⎣⎩⎭⎤−⎥−⎦∫ (2.92)

29

The total integral in equation (2.60) is divided into components along the same lines as these individual integrals:

30

j

1,2,3,4,5,jjjjjeeeeee±±±±±±=++++. (2.93)

where

()()()()()()()()24671,12,13,14,5,1sgnImsgnImsgnImsgnIm.jjljjljjljjjleFeFeFeeFQdQωωωω±±±±±±±±±±±±±

Q dQ

Q dQ

Q dQ

⎧⎫⎪⎪=−⎨⎬⎪⎪⎩⎭⎧⎫⎪⎪=−⎨⎬⎪⎪⎩ ⎭

⎧⎫⎪⎪=−⎨⎬⎪⎪⎩⎭⎧⎫⎪⎪+=−⎨⎬⎪⎪⎩⎭∫∫∫∫ (2.94)

The integrals solved above are then substituted into these equations producing

()()()()2021,10120211202,3,1secsgncossecsecsin(1)21;022sinsgncos;2sinsecsin(1)2;0221;020;2;02sjjjjjjkeknkBnjnBnBkjBejjeπθωθωπθθππθπωθωππθθθθπθπθ±±±±±±±−=+−⎧⎧=<<−⎪⎪⎛⎞⎪⎪⎨⎜⎟⎪⎝⎠⎪+−=−<<⎪⎪⎪⎩=⎨⎧⎪=−<<⎪⎪⎪⎨⎪⎪⎪=<<⎪⎪⎩⎩−=()()()()()()()()111220014,12015,120gn11;tan2sin10;tan1secsecsgnsecsin(1)2(1)sgn2sinsin2sinsgnsecsin(1)jjjjjjjxamyBxamyBkfkenkBnnfBBekωπθθθθθωωπθθπωπθθθωθθ±−−±±±±±±⎧−−+−⎛⎞⎪=⎜⎟⎪⎝⎠⎨−+−⎛⎞⎪≠⎜⎟⎪⎝⎠⎩−=−+−⎛⎞⎛⎞−−⎜⎟⎜⎟⎝⎠⎝⎠=−+−mm.2nBπ (2.95)

At this point, the wave number integral of the double integral component in the leading edge contribution (2.41) has been solved using the complex integration method. The wave angle integral in both leading edge components still must be solved.

2.3.3. Wave Angle Integral

Each of the resulting functions from the previous section, equation (2.95), is integrated with respect to the wave angle in equation (2.41). The second component in this equation is also integrated with respect to the wave angle. To simplify this 31

integration, the second component is written in terms similar to the first component (after the complex integration). The integrand of the second term (2.39) is

32

j e

()1,6,6,1njjJe+−=−− (2.96)

with

()()()2016,201cossecsecsin12jjjkenkBωθπθθ±±−−=+−. (2.97)

The total leading edge contribution (2.47) is

()()()2,,,1,11,2220202,031,11,22,sec4secsec.4nmnmLnmAkxyIkkAkJJdππππ

I dkd ζθθπθθθπ∞−−=+−++∫∫∫ (2.98)

Using the complex integration method (2.48) the wave number integration is simplified to

()()()22,,03,,121,11,2222,secsec44nmnmnmLAAkxyEEdJππππ

J d ζθθθππ−−=+++∫∫ θ . (2.99)

This leads to a natural splitting of the leading edge contribution (2.99) into two parts

,,12nmLTTζ=+ (2.100)

where

22,,031,222secsec44nmnmjjAAkTEdππππ

j J d θθππ−−=+∫∫ θ θ . (2.101)

These integrals are expressed in terms of the functions defined above (2.55), (2.93), and (2.96) as

()()()()221,2,3,4,5,221,2,3,4,5,22331,6,6,22secsec1secsec1.njjjjjjjjjnjjjEdeeeeeeeeeedJdeedππππππππθθθθθθθθ+++++−−−−−−−+−−−⎡=−++++⎣⎤−++++⎦⎡⎤=−−⎣⎦∫∫∫∫

j j

j t

(2.102)

Each part is made up of a positive and a negative contribution:

()1njjTt+−=−− (2.103)

where

2,,,1,2,3,2222,,,034,5,6,22secsecsec444secsecsec.444nmnmnmjjjnmnmnmjjjAAAteeeAAAkeeedππθθθπππ j

d θθθπππ±±±±−±±±⎛=++⎜⎝⎞+++⎟⎠∫

θ θ

(2.104)

To simplify the integration, each of these integrals is written separately using equations (2.95) and (2.97). The first integral is

()()( )

21,1,23220011,1220sec1seccossecsgn.secsin(1)2jjjjjRedkkRdnkBππππθθπθθωωθπθθ±±−±±±−=−=+−∫∫ (2.105)

33

The second integral has two forms depending on the value of j, and is

()()22,2,22122,112002102,21220secseccos2sin2sinsgn;1secsin(1)2seccos2sin2sinsgn;2.secsin(1)2jjjjRednnBBRnkBnnBBRdnk

d j

j

Bππππθθπωπθθθωπθθπωπθθθωθπθθ±±−±±±±±±−=⎛⎞⎜⎟⎝⎠==+−⎛⎞⎜⎟⎝⎠==+−∫∫∫

θ (2.106)

The third integral’s integrand only has a value at one wave angle and so is zero:

23,3,23,sec0.jjjRedRππθθ±±−±==∫ (2.107)

The fourth integral is

()()( )

24,4,23220014,1220sec1secsecsgn.secsin(1)2jjjjjRedkfkRdnkBππππθθθθωωθπθθ±±−±±±−=−=−+−∫∫ (2.108)

The fifth integral is slightly more complicated. It is

()25,5,2125,1220secsec(1)sgn2sinsin2sinsgn.secsin(1)2jjjjjRednnfBBRdnkBππππθθπωπθθθθωθπθθ±±−±±±−=⎛⎞⎛⎞−−⎜⎟⎜⎟⎝⎠⎝⎠=−+−∫∫ (2.109)

34

One of the sign functions depends solely on the value of j and whether the wave angle is positive or negative. To simplify the equation it is split into two parts; one for positive wave angles and the other for negative wave angles, i.e.

35

bj

5,5,,5,,jajRRR±±±=+ (2.110)

with

()()125,,1200105,,1220sec(1)2sin2sinsgnsecsin(1)2sec(1)2sin2sinsgn.secsin(1)2jajjjbjjnnfBbRdnkBnnfBBRdnkBπππωπθθθωθπθθπωπθθθωθπθθ±±±±±±−⎛⎞−⎜⎟⎝⎠=+−⎛⎞−⎜⎟⎝⎠=−+−∫∫ (2.111)

The last integral is

()(()

)

236,6,2322016,220sec1seccossec.secsin12jjjjjRedkRdnkBππππθθθωθθπθθ±±−±±−=−−=+−∫∫ (2.112)

Similar terms are grouped to match the equations in the Cheng and Wellicome paper. The first of these groups is

,0,,0,11,6,2444nmnmnmjjAkAAkefRRπππ j

±±=+ ± (2.113)

which becomes

()()()32201,11220seccossec11sgnsecsin(1)2jjjkefdnkBππθθωωθπθθ±±±−⎡⎤=−−−⎣⎦+−∫. (2.114)

The next three terms are simply equations from above,

36

j

j

(2.115) ,14,,25,,,35,,jjjajbelRelRelR±±±±±±===−

which become

()()()()()322101,10220112,220011,320sgnsecsec1secsin(1)2secsgnsin2sin(1)2secsin(1)2secsgnsin2sin(1)2secsin(jjjjjjjjfkelkdnkBnfBneldnBkBnfBnelBkπππωθθωθπθθπωθωθθπθπθθπωθωθθπθθ±±±−±±±±±±=−−+−⎛⎞⎜⎟⎝⎠=−+−⎛⎞⎜⎟⎝⎠=−+∫∫02.1)2jdnBπθπ−−∫ (2.116)

The positive and negative components of the leading edge contribution (2.104) are also written in terms of the Cheng and Wellicome functions:

,,,,0,,1,2,3,12,222244444nmnmnmnmnmjjjjjAAAAkAtelelelefπππππ±±±±±=+−++ j R± . (2.117)

This leads to the leading edge contributions (2.103) of

()()()()(),,,1,1,2,2,3,322,0,,1,12,2,21114411.44nnnnmnmjjjjjjnnnmnmjjjjAATelelelelelAkAefefRRππ j el

ππ+−+−+−+−+−⎡⎤=−−+−−−−+⎢⎥⎣⎦⎡⎤⎡⎤+−−+−−⎣⎦⎣⎦ (2.118)

To get the total leading edge contribution (2.100) this function is summed with j set to one and two. Based on this, the total last term of the leading edge contribution is

()()()2,,2,2,2,12,12,22,222111144nnnnmnmjjjAARRRRRRππ+−+−+−=⎡⎤⎡−−=−−+−−⎣⎦⎣Σ ⎤⎦

. (2.119)

Using the functions in equation (2.106) this equation is

()()()()()()2,2,2,21112,1120101114seccoscos2sinsin2sin1sgnsgn8secsincos22sincos2sin1sgncosnnmjjjnnmnARRnnAnBBdnBnkBBnBnπππωπωθθθ θ

ωωθππωθθθπωθωπ+−=+−+−−++⎡⎤−−=⎢⎥⎣⎦⎛⎡⎤⎛⎞⎛⎞⎜⎢⎥⎜⎟⎜⎟⎝⎠⎝⎠⎜⎢⎥−−⎜⎢⎥⎛⎞−⎜⎢⎥⎜⎟⎜⎝⎠⎣⎦⎝⎛⎞⎜⎟⎝⎠+−Σ∫()1012210seccossin2sinsgn.secsin22sinnBdnkBBππωθθθωθπωθθθ−−−−⎡⎤⎛⎞⎟⎢⎥⎜⎟⎝⎠⎟⎢⎥−⎟⎢⎥⎛⎞+⎟⎢⎥⎜⎟⎝⎠⎣⎦∫

⎞

⎠

(2.120)

To simplify the integrand, the term

11cos2sincos2sinnBnBπωθπωθ+−⎛⎞⎜⎟⎝⎛⎞⎜⎟⎝⎠

⎠ (2.121)

is studied. From the definition given in equation (2.39) and trigonometric identities, this term is

()()()11coscos1cos2sin2coscossin1sin2sin2sin2coscos1coscos2sin22sinsin2nnxamybBnnnxamyBBBnnnxamybBBnxaBπθπθπω

2

2

2

n

n

n

π

πθππθθπθππωθθπ+−⎛⎞+−+⎡⎤⎜⎟⎣⎦⎝⎠⎛⎞⎛⎞++−+⎡⎤⎜⎟⎜⎟⎣⎦⎝⎠⎝⎠⎝=⎛⎞⎛⎞⎛+−+⎡⎤⎜⎟⎣⎦⎜⎟⎝⎠⎝⎠−+()

π

⎛ ⎞

⎜ ⎟

⎝ ⎠

⎛ ⎞

⎜ ⎟⎠

⎞

⎜ ⎟

⎝ ⎠

cos1sisin22nnmyB

n θππθ⎛⎞−+⎡⎤⎜⎟⎣⎦⎝⎠

⎛ ⎞

⎜ ⎟

⎝ ⎠

. (2.122)

37

If n is odd, the cosine terms are zero, leaving

()()11coscossin1sin2sin2sin2cossin1sincos2sin22sinnnnxamyBBBnnnxamyBBBπω

2

2

n

n

πθππθθπθπππωθθ+−⎛⎞⎛⎞+−+⎡⎤⎜⎟⎜⎟⎣⎦⎝⎠⎝⎠⎝=⎛⎞⎛⎞⎛−+−+⎡⎤⎜⎟⎣⎦⎜⎟⎝⎠⎝⎠

⎛ ⎞

⎜ ⎟⎠

⎞

⎜ ⎟

⎝ ⎠

(2.123)

which simplifies to

11cos2sin1cos2sinnBnBπωθπωθ+−⎛⎞⎜⎟⎝⎠=−⎛⎞⎜⎟⎝⎠. (2.124)

On the other hand, if n is even, the sine terms are zero, leaving

()()11coscoscos1cos2sin2sin2coscos1coscos2sin22sinnnnxamyBBBnnnxamyBBBπω

2

2

n

n

πθππθθπθπππωθθ+−⎛⎞⎛⎞+−+⎡⎤⎜⎟⎜⎟⎣⎦⎝⎠⎝⎠⎝=⎛⎞⎛⎞⎛+−+⎡⎤⎜⎟⎣⎦⎜⎟⎝⎠⎝⎠

⎛ ⎞

⎜ ⎟⎠

⎞

⎜ ⎟

⎝ ⎠

(2.125)

which simplifies to

11cos2sin1cos2sinnBnBπωθπωθ+−⎛⎞⎜⎟⎝⎠=⎛⎞⎜⎟⎝⎠. (2.126)

Therefore, overall this function is equivalent to

()11cos2sin1cos2sinnnBnBπωθπωθ+−⎛⎞⎜⎟⎝⎠=−⎛⎞⎜⎟⎝⎠. (2.127) 38

So the equation (2.120) simplifies to

()()()()()2,2,2,2112,11200101122014seccossin2sinsgnsgn8secsin2seccossin2sinsgnsgn.secsin2nnmjjjnmARRnAnBdnBkBnBdnkBππππωθθθωωπθθπωθθθωωθπθθ+−=−+−−+−−⎡⎤−−=⎢⎥⎣⎦⎛⎛⎞⎜⎜⎟⎝⎠⎜⎡⎤−⎣⎦⎜−⎜⎝⎞⎛⎞⎟⎜⎟⎝⎠⎟⎡⎤+−⎣⎦⎟+⎟⎠Σ∫∫

θ (2.128)

Writing this in terms of Cheng and Wellicome functions produces

()(2,,2,2,1,22,221148nnmnmjjjAAn ) RRefBπ+−=⎡⎤−−=+⎢⎥⎣⎦Σ ef (2.129)

where

()()()()()()121,211200102,211220seccossin2sinsgnsgn;1secsin12seccossin2sinsgnsgn;2.secsin12jjnBefdjnkBnBefdjnkBπππωθθθωωθπθθπωθθθωωθπθθ−+−−+−−⎛⎞⎜⎟⎝⎠⎡⎤=−⎣⎦+−⎛⎞⎜⎟⎝⎠⎡⎤=−⎣⎦+−∫∫

=

=

(2.130)

39

To summarize, the total leading edge contribution is

()()()()()()()(),,,1,11,11,21,21,31,322,12,12,22,22,32,3,01,11,12,12,1,1,22,211141111148nnnnmnmLnnnnnnmnmAelelelelelelelelelelelelAkefefefefAnefefBζππ+−+−++−+−+−+−+−⎡=−−+−−−−+⎣⎤+−−+−−−−+⎦⎡⎤+−−+−−⎣⎦⎡⎤++⎣⎦

−

(2.131)

using the following integrals

()()()()()322101,10220112,220011,320sgnsecsec1secsin(1)2secsgnsin2sin(1)2secsin(1)2secsgnsin2sin(1)2secsin(jjjjjjjjfkelkdnkBnfBneldnBkBnfBnelBkπππωθθωθπθθπωθωθθπθπθθπωθωθθπθθ±±±−±±±±±±=−−+−⎛⎞⎜⎟⎝⎠=−+−⎛⎞⎜⎟⎝⎠=−+∫∫()()()()()()()()0232201,11220121,2112002,2111)2seccossec11sgnsecsin(1)2seccossin2sinsgnsgn;1secsin12seccossinsgnsgnjjjjjdnBkefdnkBnBefdjnkBnefππππθπθθωωθπθθπωθθθωωθπθθθθωω−±±±−−+−+−−⎡⎤=−−−⎣⎦+−⎛⎞⎜⎟⎝⎠⎡⎤=−⎣⎦+−⎡⎤=−⎣⎦∫∫∫ =

()102202sin;2secsin12jBdjnkBππωθθπθθ−−⎛⎞⎜⎟⎝⎠=+−∫ .

(2.132)

These integrals do not have an analytical solution; they must be numerically integrated. Several of these integrals contain the even auxiliary function from equation

40

(2.91). This function also does not have an analytical solution and either a series expansion or numerical integration must be employed.

2.3.4. Evaluation Grid

As was mentioned in equations (2.45) and (2.46), the wave height is given by

()()()()(),,,,111,,sin.2NMnmLnmLnmNnngxygxyxaynAxyBBρζρζζπ===

, ⎡⎤=− − ⎣⎦⎡⎤−+⎢⎥⎣⎦ΣΣΣ (2.133)

In the previous section it was determined that this equation cannot be solved analytically; therefore, a numerical solution is employed. To approximate this function numerically, a grid of solution points is created. Two indices are already in use; m corresponds to the longitudinal pressure strip in question and n corresponds to the frequency of the sine wave in the pressure. To define where on the surface the height is being solved for, two more indices are used: i corresponds to the longitudinal (x) direction and j corresponds to the transverse (y) direction. Using these indices the x and y coordinates become

121.2ijLxiIByjJ⎛⎞=−⎜⎟⎝⎠⎛⎞=−⎜⎟⎝⎠ (2.134)

I is the number of longitudinal points in the wetted area and J is the number of transverse points.

41

Fig. 7 Sinusoidal Pressure Hull Discretization

i = 1

x = x1

i = 2

x = x2

x = -L

x = 0

x = -a

m = 1

m = 2

m = M

j = J

y = yJ

j = 2

y = y2

y = B

j = 1

y = y1

y = 0

i = I

x = xI

y,η

x,ξ

a

Similarly, the pressure coefficient from the last term of the wave height (2.133) is written in terms of the longitudinal index:

(),ninAxA= i . (2.135)

At this point a tensor containing the wave height from the leading edge contribution of each pressure strip and pressure frequency for each of the grid points is created:

(),,,,,,,nmLijnmijnmgxyEAρζ=. (2.136)

It has been noted that the trailing edge contribution is expressed by shifting the leading edge equation by a patch length (2.43), i,e,

12ixaaia⎛⎞+=−+⎜⎟⎝⎠. (2.137)

This is equivalent to shifting the index of the longitudinal coordinate:

42

43

i

1ixax−+=. (2.138)

The same logic is used to express each of the subsequent pressure strips as the first strip offset:

()()()111121.iiimxamaiamxamx+−⎛⎞+−=−+−⎜⎟⎝⎠+−= (2.139)

These simplifications are then substituted back into the wave height equation (2.133):

()(,,,1,1,,1,,,111sin2NMNijnmnimjnimjnijnmnngAEEAyBπρζ+−−===

B) ⎡⎤=−− + ⎢⎥⎣⎦ΣΣΣ. (2.140)

This equation is further simplified by combining the n summations. Additionally, the case where the transverse position is outside the half breadth of the craft (2.26) is now added back in, resulting in

()()(),,1,1,,1,,11,,,,1,1,,1,,11;sin2;NMnmnimjnimjnmjnijijNMjnmnimjnimjnmjAEEByBnAyBgByBAEEByπρζ+−−==+−−==⎧⎛−⎪⎜⎝⎪−≤≤⎪⎞⎡⎤⎪−+=⎨⎢⎥⎣⎦⎠⎪⎪≤−⎪−≤⎪⎩ΣΣΣΣ

⎟ . (2.141)

The m sum is combined with the last term by including a condition when m = i,

()()()(),,1,1,,1,,11,,,1,1,,1,0,,,,1,1,,1,,11;sin2;NMnmnimjnimjnmmijijninjnjnijNMjnmnimjnimjnmjAEEByBngAEEAyBByBAEEByπρζ+−−==≠+−−==⎧⎡⎪⎢−⎪⎢⎪⎢−≤≤⎣⎪⎪⎤⎛⎞=⎡⎤⎨+−−+⎥⎜⎟⎢⎥⎪⎣⎦⎝⎠⎦⎪≤−⎪−⎪≤⎪⎩ΣΣΣΣ. (2.142)

According to the discretized longitudinal position equation (2.134), the term En,1,0,j corresponds to a location ahead of the craft. Therefore, the wave height due this term is obviously zero. A new wave height basis function tensor is generated that incorporates each of these changes.

(),1,1,,1,,*,,,,1,1,; or or sin; and -2nimjnimjjjnmijnjjjEEimyByEnEyBimByBπ+−− B

B

−≠<−⎧⎪=⎨⎡⎤−+=≤≤⎪⎢⎥⎣⎦⎩

>

, , ,

. (2.143)

Substituting back into the total wave height (2.142) produces

. (2.144) *,,11NMijnmnmijnmgAEρζ===ΣΣ

This equation is linear in the pressure coefficients. It can be written in matrix form by introducing two new indices: p corresponds to the pressure strip and frequency and k corresponds to the location on the surface, i.e.

()()11.pmNkiJj

n =−+=−+ (2.145) 44

The wave height equation then becomes

45

p k

(2.146) *,1MNkppgAEρζ==Σ

which is written in matrix form as

*gρζ= A E. (2.147)

In summary, the integrals in equation (2.132) are numerically integrated over a grid of points for the leading edge of the first pressure strip. The grid is then shifted to generate the entire basis function tensor. Using the basis function tensor and pressure coefficients the wave height is determined.

2.4. CONSTANT PRESSURE PATCHES

A second method to define the pressure patches is to assume a constant pressure over the entire patch [Scullen]. Multiple pressure patches with varying pressures are then created over the wetted area of the craft. The wave height is linear in each of the contributing heights, so the heights from the individual pressure patches are added together to generate the total wave height, i.e.

()(,11,NMnmnm

, ) xyζζ===ΣΣ x y . (2.148)

Scullen and Tuck present the following wave height

()(),,,,,nmnmnmzSP

,0 xyGxydg

d ζξηξρ=−−∫∫ η (2.149)

created by such a pressure patch, where

()()2cossin202020,,Resec2seikxykzkeGxyzdkdkkπθθπ

c2

θθπθ∞−++−=−∫∫. (2.150)

The partial derivative of this equation with respect to height, evaluated on the undisturbed free surface (z = 0), is

()()2cossin2022020,,0Resec2sikxyzkkeGxydkdkkπθθπ ec

θθπθ∞−+−=−∫∫. (2.151)

Therefore, the wave height due to a single constant pressure patch is

()()(),2,02,22020cossin,Resec2.nmnmnmikxySPkkxygkkedddkdππξθηθζθ

sec ρπθξηθ∞−−−+−⎡⎤⎣⎦=−∫∫∫∫g (2.152)

Fig. 8 Constant Pressure Hull Discretization

x = -L

x = 0

x = -a

m = 1

m = 2

m = M

b

y,η

y = 0

y = B

i = 1

x = x1

i = 2

x = x2

i = I

x = xI

j = 1

y = y1

j = 2

y = y2

j = J

y = yJ

n = N

n = 2

n = 1

y = b

x,ξ

a

As with the sinusoidal pressure method, a double integral over the wetted area of the craft must be evaluated to determine the surface height.

46

2.4.1. Surface Integral

From Fig. 8 the surface integration is

()()()()()(),1cossincossin1nmambnikxyikxySambneddeξθηθξθηθ d d ξηξ−−−−+−−−+−⎡⎤⎡⎤⎣⎦⎣⎦−−=∫∫∫∫ η . (2.153)

To evaluate this integral two variable substitutions are made:

**.xyξξηη=−=− (2.154)

The surface integral evaluates to

()()()()()()()()()(),cossincossincossin22cossincossin22sincossincos.sincossincosnmikxySikxamaybnikxamaybnbikxamybnikxamybnbeddeekkeekkξθηθθθθθθθθξηθθθθ

θ

θθθθ−−+−⎡⎤⎣⎦−+−+−−+−+−+⎡⎤⎡⎣⎦⎣−++−−++−+⎡⎤⎡⎣⎦⎣=−++−∫∫

⎤⎦ ⎤⎦

(2.155)

The single patch wave height (2.152) becomes

()()()()()()()()()()22,0,22200cossincossincossincossinsec1,Re2sincossec.nmnmikxamaybnikxamaybnbikxamybnikxamybnbPkxygkkkeeeeππ

dkd

θθθθθθθθζρπθθ

θ

θθ∞−−+−+−−+−+−+⎡⎤⎡⎣⎦⎣−++−−++−+⎡⎤⎡⎤⎣⎦⎣⎦=−−⎡−⎣⎤−+⎦∫∫g ⎤⎦

(2.156)

The surface integration is accomplished. Two integrations remain in this method, the wave angle and wave number integrations. The first one to evaluate is the wave number integral.

47

2.4.2. Wave Number Integral

To simplify the single patch wave height, a new function is created based on the four terms of the surface integral. Essentially, a height component must be found for each of the four corners of the pressure patch. The wave height equation becomes

()()(()(),,0000,,,,,.nmnmP ) xyxamaybnxamaybnbgxamybnxamybnbζζζρζζ=+−−−+−−+⎡⎣−+−++−+⎤⎦ (2.157)

where

()()()22cossin0022200sec11,Re2sincossecikxykxyekkkπθθπθ

dkd ζθπθθθ∞−+−=−−∫∫. (2.158)

Evaluating the wave number integral produces

()()()2022glog2Hsin1,2sincosTTTTxydπππζθπθθ−+−=−∫ (2.159)

with

200secsecsinTkxkyθθ=+ θ , (2.160)

g the even auxiliary function defined in equation (2.91), and H the Heaviside step function, defined as

()0;0H1;0zzz<⎧=⎨>⎩. (2.161)

2.4.3. Wave Angle Integral

As with the sinusoidal pressure patches, the wave angle integral must be evaluated numerically. To simplify this integration the corner integral (2.159) is divided into two parts, a local and a far field.

48

()()()0,,LF , xyxyxζζζ=+ y (2.162)

The local field corresponds to the area directly under the craft and does not contribute to the wake created by the craft; it does have an effect near the craft, however. This term is

()()()()222glog2HHsin1,2sincosLTTTxTxydπππζθπθθ−+−−⎡⎤⎣⎦=−∫. (2.163)

The far field on the other hand contains the craft wake, but the integrand is much simpler:

()()22Hsin,sincosFxTxyππ

d ζθπθθ−=∫. (2.164)

The total wave height is found by summing the far field component at each corner of the pressure patch and accounting for the local field contribution which will be discussed in section 2.5.2.

2.4.4. Evaluation Grid

As with the sinusoidal pressure patches, the wave angle integral is solved numerically. Therefore, a grid over which to solve the integral is created. The same grid as the sinusoidal pressure method is used with one difference. Previously, n corresponded to the frequency of the pressure; now it corresponds to the transverse pressure patch. The location is discretized in the same way as in the previous method (2.134). Similar to the sinusoidal pressure method, a basis tensor consisting of the contributions of each pressure point at each index location is found:

()()()(),,,0000,,,,nmijijijijijExamaybnxamaybnbxamybnxamybnbζζζζ=+−−−+−−+−+−++−+ .

(2.165)

This method requires a wave height contribution from each of the corners of the pressure patch. Each of these contributions is identical except for the offset. This offset is

49

()()121211.iijjxamaiamxamaaiamaybnbjbnybnbbjbnb⎛⎞+=−+⎜⎟⎝⎠⎛⎞+−=−+−⎜⎟⎝⎠−=−−−+=−−+ (2.166)

These can also be written in terms of the index as

11.iimijjnjjxamxxamaxybnyybnby−

i m

n

−+−−++=+−=−=−+= (2.167)

The basis tensor (2.165) then becomes

()()()(,,,01011001,,,,nmijimjnimjnimjnimjnExyxyxyxyζζζζ−+−−+−+−−−−+=−−+ ). (2.168)

As with the previous method, only the first pressure patch contribution needs to be calculated. The other patches are identical but offset. Substituting this basis tensor back into the single pressure patch wave height equation (2.157) produces

(),,,nmijnmnmijgxyPEρζ= , , ,

, , ,

. (2.169)

The total wave height (2.148) therefore becomes

. (2.170) ,,11NMijnmnmijnmgPEρζ===ΣΣ

The system is linear in the pressure coefficients and is written in matrix form similar the sinusoidal pressure method, equations (2.145), (2.146), and (2.147).

In this method equation (2.164) must be numerically integrated. The integration is performed over a grid of points to produce the contribution of one corner of a pressure patch. The total patch contribution is found by shifting the first corner to the other three

50

51

corners and summing the results. This first patch is then shifted to produce the entire basis function tensor. The pressure amplitudes together with this tensor determine the total wave height.

2.5. NUMERICAL ISSUES

Both the sinusoidal and the constant pressure methods involve a numerical integration of the wave angle integral. These integrals all contain features such as singularities and rapid oscillations which make the numerical integration difficult to perform.

2.5.1. Sinusoidal Pressure Method

The integrals for the sinusoidal pressure method (2.132) are

()()()()()3220101,1220112,220011,320sgnsecsec1secsin(1)2secsgn2sin2sin(1)secsin(1)2secsgn2sin2sin(1)secsin(jjjjjjjjkfkeldnkBnnfBBeldnkBnnfBBelkπππωθθωθπθθπωπθωθθθπθθπωπθωθθθθ±±±−±±±±±±=−−+−⎛⎞⎜⎟⎝⎠=−+−⎛⎞⎜⎟⎝⎠=−+∫∫()()()()()()()()02322101,122011121,2200112,21)21sgnseccossec1secsin(1)2secsgnsgncossin2sin;1secsin12secsgnsgncossinjjjjjdnBkefdnkBnBefdjnkBnefππππθπωθθωθπθθπωθωωθθθπθθθωωθ−±±±−−+−+−−⎡⎤−⎣⎦=−−+−⎛⎞⎡⎤−⎜⎟⎣⎦⎝⎠==+−⎡⎤−⎣⎦=∫∫∫()102202sin;2secsin12jBdjnkBππωθθπθθ−−⎛⎞⎜⎟⎝⎠=+−∫ .

(2.171)

It should be noted that the denominator is identical in each of these integrands:

20secsin(1)2jjnDkBπθθ=+−. (2.172) 52

The integrands also contain only a few common numerators:

()()()()()()()32010111321011111sgnsecsecsec2sgn2sin2sin11sgnseccossecsec2sgnsgncos.sin2sinLNkfknnLNfBBFNknFNBωθθωπωπθωθθωθθωπωθωωθθ±±±±±±±±±−+−=⎛⎞=⎜⎟⎝⎠⎡⎤=−⎣⎦⎛⎞⎡⎤=−⎜⎟⎣⎦⎝⎠ (2.173)

To cut down on the number of calculations performed during the numerical integration, the leading edge wave height equation (2.131) is divided into three terms:

()()(,,0,,,2448nmnmnmnmLAAkAn) IIIIIIBζππ=++. (2.174)

The first term in this equation is

()()()()()()1,11,11,21,21,31,32,12,12,22,22,32,3111111nnnnnnIelelelelelelelelelelelel+−+−+−+−+−+−=−−+−−−−++−−+−−−−+ .

(2.175)

This is written in terms of the numerators and denominator above as

()()()()()()()()01222102111122111122111122111122.nnnnnnnnLNLNLNLNIDLNLNLNLNdDLNLNLNLNDLNLNLNLNdDππθθ+−+−−+−+−+−+−+−+−⎡−−+−−=⎢⎢⎣⎤−−+−−++⎥⎥⎦⎡−−−−++⎢⎢⎣⎤−−++−−+⎥⎥⎦∫∫ (2.176)

53

The second term is

()()1,11,12,12,111nnIIefefefef +−+=−−+−− −; (2.177)

or in terms of the common numerators and denominators

()()2122111111nnFNFNFNFNIIdDDππθ+−+−−⎡⎤−−−−+=+⎢⎥⎢⎥⎣⎦∫. (2.178)

The third term is

1,22,2IIIefef=+; (2.179)

or in the common numerators and denominators

20120222FNFNIIIddDDππθθ−=+∫∫. (2.180)

Combining the integrals in this way cuts down on the number of times the integration routine calculates common terms. It also allows the integration routine to only be more accurate where needed and reduce the number of integrating points thereby speeding up the calculation. Typical combined integrands are shown below.

54

-pi/2-pi/40pi/4pi/2-30-20-100102030Iθ2θ−4θ−04θ2θ

-pi/2-pi/40pi/4pi/2-30-20-100102030IIθ2θ− 4θ−0 4θ2θ -pi/2-pi/40pi/4pi/2-30-20-100102030IIIθ2θ−4θ−04θ 2θ

Fig. 9 Sinusoidal Pressure Combined Integrands

k0 = 3.08, x = 0.25, y = 0.14, m = 1, n = 1, a = 0.1, B = 0.5

There are still a few issues preventing accurate numerical integration of these functions. The first is the rapid oscillation in II as it approaches the integration limits. To aid in the integration of this oscillation a variable substitution is made:

tantθ=. (2.181)

This substitution applied to the II integral (2.178) leads to

()()****12111111nnFNFNFNFNIIdDD+−+−∞−∞⎡⎤−−−−+=+⎢⎥⎢⎥⎣⎦∫ t (2.182)

55

with the numerator and denominator being

()()()()*210201211sgn1cos11(1)21.1jjFNtktnDkttBxamyBtt

2

1 ωωπω±±±⎡⎤⎡=−++⎣⎦⎣=++−+−+=+m

±⎤⎦

(2.183)

This new integrand is shown below. -505-505IIt

Fig. 10 Sinusoidal Pressure II Integrand with Slow Oscillations

k0 = 3.08, x = 0.25, y = 0.14, m = 1, n = 1, a = 0.1, B = 0.5

The oscillations are much more reasonable. There is a new problem, however, the integration limits are now positive and negative infinity. The integrand rapidly approaches zero as t approaches positive or negative infinity. The integration is safely truncated with minimal inaccuracies in the result.

Another problem exists; the integrand of each of the functions becomes infinite at some points. To find where these singularities are, the value of the wave angle that sets the denominator equal to zero needs to be found:

2**0secsin(1)0.2jjnDkBπθθ=+− = (2.184)

56

Instead of solving this equation directly, it is easier to solve it in terms of the transformed variable used to stretch the II integrand (2.181):

**201(1)2jjnDkttB

0

π=++−=. (2.185)

The solution where this singularity occurs is

2*011122ntBkπ⎛⎞=±−++⎜⎟⎝⎠. (2.186)

The wave angle of the singularity is therefore

2*1011tan122nBkπθ−⎡⎤⎛⎞⎢⎥=±−++⎜⎟⎢⎥⎝⎠⎢⎥⎣⎦. (2.187)

By inspection, when j = 1 in the denominator the sign of the wave angle must be positive and when j = 2 the sign must be negative to make the total denominator zero.

The location of the singularity has been found, now the singularity must be dealt with. To remove the singularity from the integration, the integrands are modified. In general,

()22****11*2**1lnNNNNddDDDDθθθθθθθθθθθθθθθθθθθθθθ====⎡⎤⎛⎞−⎢⎥=−+⎜⎜−−⎢⎥⎝⎠⎣⎦∫∫ ⎟⎟

. (2.188)

The second term in the integrand effectively cancels out the singularity of the first term. To not add any net change to the equation, the integral of the second term is added back to the total equation. This technique is used to remove the singularity from each of the integrals in both the original (wave angle) and transformed variable. To apply this technique the derivative of the denominator must be computed. In the original variable the derivative is

57

()2302cossecDkθθθ=−; (2.189)

in the transformed variable the derivative is

()202121ktDtθ+=+. (2.190)

The modified integrands with the singularities removed are shown below.

-pi/2-pi/40pi/4pi/2-10-505Iθ2θ−4θ−04θ2θ -10-50510-4-3-2-101IIt -pi/2-pi/40pi/4pi/2-20246IIIθ2θ−4θ−0 4θ2θ

Fig. 11 Sinusoidal Pressure Final Integrands

k0 = 3.08, x = 0.25, y = 0.14, m = 1, n = 1, a = 0.1, B = 0.5

The integrands still contain discontinuities, but they are now either due to different integration ranges, i.e. I and III contain two integrals one for negative wave angles and the other for positive wave angles, or to the sign function contained in the integrands. 58

The first discontinuity type is not an issue because over each integration range the function is continuous. The sign discontinuity could be removed by splitting the integration into two separate integrals at this point. This introduces several coding issues though; locating the discontinuity and determining which segment the singularity is in to name a couple. Because of these issues the sign discontinuity is resolved by simply using many integration points. While this is not the most elegant method, it does not introduce significant integration errors.

2.5.2. Constant Pressure Method

The constant pressure method has two terms to evaluate, equations (2.163) and (2.164); they are

()()()()()()22222glog2HHsin1,2sincosHsin,sincosLFTTTxTxydxTxydπππππζθπθθζθπθθ−−+−−⎡⎤⎣⎦=−=∫∫ (2.191)

with

200secsecsinTkxkyθθ=+ θ . (2.192)

The local field integral is treated specially. As Scullen and Tuck show, the local field contribution to the surface height is mainly directly under the hull. It appears as a discontinuity in the far field contribution that is constant in the longitudinal direction. It also contains the even auxiliary function (2.91) which must be computed either from an interpolation or a series expansion. Therefore, instead of taking the computational time and effort, the local field contribution is approximated as the value of the far field at the leading edge of the wetted area. The total wave height at this location is effectively zero. 59

Therefore, the local field must cancel the far field out. The local field is assumed constant in the longitudinal direction for the length of the wetted area.

All of the actual wave height in the constant pressure method comes from the far field contribution. This integral must still be computed. The far field function is much simpler than the previous integrals. It does not have any complex functions in the integrand and only contains a single integral to evaluate. It does have the same issues as the sinusoidal pressure integrals. The integrand oscillates as the wave angle approaches the integration limits and it also contains a singularity. These problems are solved using the same methods developed for the sinusoidal pressure method.

-pi/2-pi/40pi/4pi/2-2002040ζFθ2θ−4θ−04θ2θ

Fig. 12 Constant Pressure Original Far Field Integrand

k0 = 3.08, x = 0.25, y = 0.14

The integrand oscillation is countered by a variable transformation (2.181). The far field integral then becomes

()()Hsin,FxTxytζπ∞−∞=∫ dt (2.193)

with

60

(201Tktxyt=++ ). (2.194)

This integrand is shown below. -10-50510-10-50510ζFt

Fig. 13 Constant Pressure Far Field Integrand with Slow Oscillations

k0 = 3.08, x = 0.25, y = 0.14

This integrand still has a singularity at t = 0. The singularity is removed using the method from the sinusoidal pressure patches (2.188). Using this method the far field integral becomes

()()()0Hsinsin,FxTkxxytζπ∞−∞−=∫ dt. (2.195)

The second term in the integrand is odd so its integral is zero and nothing needs to be added to the integral to cancel it out. This integrand is shown below.

61

-100-50050100-1-0.500.51ζFt

Fig. 14 Constant Pressure Far Field Final Integrand

k0 = 3.08, x = 0.25, y = 0.14

The singularity has been removed; due to the scale of the integration limits the oscillation is more pronounced. Many integration points must be used, but this function can be accurately numerically integrated.

62

3. METHOD VALIDATION

To validate the methods presented in the previous section several tests are performed. These tests compare the sinusoidal and constant pressure methods derived above and where possible include or reference outside results.

3.1. SINUSOIDAL PRESSURE BASIS FUNCTION

The first test is taken from Cheng and Wellicome’s paper presenting the sinusoidal pressure method. The test utilizes the simplest pressure function available to this method, a single sine wave with unit amplitude,

1sin2Pπ y ⎡⎤⎛=+⎜

⎞⎟⎢⎥⎝⎠⎣⎦. (3.1)

The beam and length of this “hull” are both one. The pressure is constant in the longitudinal direction, so only one pressure strip is used, as shown in figure 1.

63

00.20.40.60.81-0.500.5xy0.20.40.60.8

Fig. 15 Sinusoidal Pressure Basis Function Test Pressure

The fundamental wave number, representing craft speed, is 3.0779. The density of the water and acceleration due to gravity are both set to one. This test is shown in figure 4 of Cheng and Wellicome’s paper. Both the sinusoidal and constant pressure methods are implemented to attempt to reproduce this figure. -20246810-505xyConstant Pressure MethodSinusoidal Pressure Method-1-0.500.51

Fig. 16 Sinusoidal Pressure Basis Function Test Wake

The gross aspects of the two methods match each other quite well and also match the results from Cheng and Wellicome. The results are very similar in the “wake” region of the surface. The peaks and valleys of the oscillations occur in the same places and are approximately the same amplitudes. On the other hand, there are significant differences

64

in the two methods. It is expected that the sinusoidal pressure method would perform better in this test because the input pressure is exactly one of its basis functions. However, that method has significantly more noise both ahead of the applied pressure and between the pressure and the limits of the wake. Ahead of the craft the surface height should obviously be zero because the craft’s speed carries the disturbance downstream. Similarly, the double rows of spikes spreading out at nearly 45 degrees behind the pressure disturbance occur ahead of the predicted wake and are noise created by this method. This noise does not create a problem using this method though because it is outside of the area of interest and does not corrupt the key results. The constant pressure method has no problem with noise in its results.

A much larger problem exists with the sinusoidal method directly under the applied pressure, where the “hull” is. The sinusoidal pressure method predicts that the water immediately rises as a pressure is applied to the surface as seen in Fig. 17. 00.20.40.60.81-0.500.5xyConstant Pressure MethodSinusoidal Pressure Method-1-0.500.51

Fig. 17 Sinusoidal Pressure Basis Function Test Hull

This rise is quite dramatic when viewed in the total wake surface, Fig. 16. This rise cannot happen physically. Approximately half way along the length of the pressure the water surface drops to match the surface predicted by the constant pressure method. This

65

66

result is not predicted in the Cheng and Wellicome paper and must be a numerical anomaly in the calculations, but no reason has been found for it. This result is much more significant than the noise of the sinusoidal method. The goal of this research is to match a hull shape to a predicted wake. If the hull the sinusoidal method predicts is incorrect in this simple example, then the method is not a valid model for predicting more complex hull shapes.

Another point that should be mentioned is the time to compute the results using each of these methods. The sinusoidal pressure method takes orders of magnitude longer to complete the calculations than the constant pressure method. This timing makes sense because the sinusoidal pressure method computes many more integrals than the constant pressure method. Additionally, the integrals in the sinusoidal method often require an additional integration at each point in the integrand. Even though this interior integral is accomplished with an interpolation, it still slows the integration down.

Overall, both methods predict very similar wake shapes. The sinusoidal pressure method has more numerical integration noise, but the noise is in parts of the surface that do not interfere with the hull and wake surfaces. The sinusoidal pressure method also takes significantly longer to compute. The major downfall of the sinusoidal pressure method, though, is the rise in the wave height near the forward section of the hull. This rise does not match any physical phenomenon (such as a bow wave) in either size or shape; additionally it does not match the results in Cheng and Wellicome. The rise makes this method unreliable for predicting a hull shape given a pressure definition.

3.2. SIMPLIFIED REALISTIC HULL PRESSURE

Scullen and Tuck also present a test case for a pressure distribution. This test case represents a simplified function for the pressure a planing hull creates. The pressure is given by

( ) 22141Px⎛⎞=−−⎜⎟⎝⎠

y . (3.2)

This “hull” has a length of two and a beam of one. 00.511.52-0.500.5xy1234

Fig. 18 Simplified Realistic Hull Test Pressure

For this test the fundamental wave number is one. As with the previous test gravitational acceleration and the density of the water are both set to one.

67

-20246810-505xyConstant Pressure MethodSinusoidal Pressure Method-10-505

Fig. 19 Simplified Realistic Hull Test Wake

Once again, the major aspects of the results from both the constant and sinusoidal pressure methods match. They both show a series of peaks and valleys in the approximate shape of an actual boat wake. The area inside the major wake has smaller disturbances. These results are very promising for the goal of using these methods to predict an actual wake. The sinusoidal method has several of the same issues that plagued the sinusoidal pressure basis function test. The noise ahead of the boat is not as apparent, but the double row of spikes behind the applied pressure is present as before, Fig. 16, and in this test intersects the actual wake making the issue more problematic. The rise in the surface under the applied pressure is also still present. This rise matches the characteristics of the previous test, a rise in the front half of the pressure becoming a depression further along the pressure.

68

00.511.52-0.500.5xyConstant Pressure MethodSinusoidal Pressure Method-10-505

Fig. 20 Simplified Realistic Hull Test Hull

The area directly under the applied pressure using the constant pressure method matches the results obtained by Scullen and Tuck, figure 9. As with the sinusoidal pressure basis function test, the constant pressure method completed the calculations many orders of magnitude faster than the sinusoidal pressure method.

The results of this test are essentially the same as the previous test. The constant pressure method performs quite well and the sinusoidal pressure method has several problems that make it an unreliable candidate to predict a hull and wake shape. An additional promising result is the wake behind this pressure is beginning to resemble a boat’s wake.

3.3. SINUSOIDAL PRESSURE CENTERLINES

The previous tests explored the general shape of the wakes generated using the two pressure methods. Cheng and Wellicome also present centerlines of the surfaces produced by the sinusoidal pressure basis function test run at different beam to length ratios. The pressure for these tests is the same as in the sinusoidal pressure basis function test (3.1) except the beam is varied in each test, i.e.

69

1sin2yPBπ⎡⎤⎛=+⎜

⎞⎟

⎢⎥⎝⎠⎣⎦ (3.3)

where B is the beam. The length and speed in this test are also the same as the sinusoidal pressure basis function test. Cheng and Wellicome present the results of these tests in figure 3. Their results are also presented in the Fig. 21. The most obvious feature of each of these plots is the pressure rise in the sinusoidal pressure method at the leading edge of the pressure. In each test, the surface this method predicts rises dramatically before dropping to a more reasonable value halfway along the pressure, x = -0.5. Other details are also different in the three methods. The first test, B = 0.4, is inconclusive. Neither method matches Cheng and Wellicome’s results. The basic form of the three results is similar; a depression under the “hull” followed by a rise and second depression aft of the boat. The details of the three are quite different, however. The second test, B = 1, provides more insight. Again, the basic forms of the results are similar. The amplitudes, though, of both the sinusoidal and constant pressure methods are not as high as Cheng and Wellicome’s results predict; they are approximately equal to each other, however. This test corresponds to the sinusoidal pressure basis function test presented above. In the sinusoidal pressure basis function test, the amplitude of the waves is around one. This is also true in Cheng and Wellicome’s results for this test. Their centerline, on the other hand, has an amplitude of over two. This implies the amplitude of their centerline tests may not be accurate. The other aspects of this test match Cheng and Wellicome’s results well. The location of the second peak does not exactly match the peak of Cheng and Wellicome, but it is close in each method. The last test, B = 10, has similar results. The Cheng and Wellicome amplitude is higher, but the peaks line up

70

better than the second test. Once again, the sinusoidal and constant pressure methods results match quite well. From Cheng and WellicomeSinusoidal PressureConstant Pressure -3-2-10-1012xz

B = 0.4 -3-2-10-2-101234xz

B = 1 -3-2-10-20246xz

B = 10

Fig. 21 Sinusoidal Pressure Basis Function Centerlines

Overall, the results from the methods tested do not match Cheng and Wellicome’s results exactly. It is expected that the sinusoidal pressure method would. Some difference is anticipated and acceptable in the constant pressure method. In each test, the

71

72

amplitudes do not match those predicted by Cheng and Wellicome. The sinusoidal and constant pressure method amplitudes do match each other reasonably well, however. The problem may be in the amplitude predicted by Cheng and Wellicome. Each method predicts approximately the same frequency for the waves; wave peaks are close enough to be acceptable. The surface rise in the sinusoidal pressure method, as in the previous tests, is a fatal problem. The constant pressure method is a valid method to predict hull and wake shapes given a pressure distribution.

73

4. METHOD ISSUES

In addition to the test case inaccuracies discussed in the previous chapter, there are several other issues which may complicate using each of the methods to compute a boat wake. Each of these issues deals with the different pressures the method is based around and involves symmetry and oscillation frequencies.

4.1. CONSTANT PRESSURE METHOD SYMMETRY

The first issue involves the fact that boats are generally symmetric. The symmetry assumes that the port and starboard sides of the hull are identical and also that the center of gravity is directly over the centerline of the hull. These two assumptions lead to a symmetric pressure distribution. Because the pressure is symmetric it is safe to only analyze one side of the wake response and assume the other side is identical. The constant pressure method utilizes constant pressure elements as the basis of the applied pressure. Multiple pressure basis functions (at different locations) are added together to produce the total pressure. A problem arises when only computing the response of half of the pressure patches, however. The response of each pressure patch is symmetric about the center of the pressure patch, but the pressure patch is offset from the centerline of the hull making the pressure (and the response) non-symmetric about the centerline of the hull (Fig. 22). Because this response is not symmetric it is no longer accurate to only examine one side of the wake.

0246810-6-4-20246xy

Fig. 22 Single Constant Pressure Patch Wake

One solution to this problem is to not take advantage of the symmetry and compute responses from pressure patches over the entire hull. This method must be employed if the hull is not symmetric for any reason. On the other hand, if the hull is symmetric there are corresponding pressure patches on either side of the centerline. These two patches can be analyzed to produce a combined response. 0246810-6-4-20246xy

Fig. 23 Corresponding Constant Pressure Patch Wakes

This response is now symmetric; however, the advantage in time saving from only computing responses for one half of the pressure patches has been lost. The response needs to be calculated for only one half of the wake, but for each of the pressure patches.

74

Another solution is to only use the pressure patch on one side, and reflect the wake created across the centerline back into the side being calculated. This reflection actually comes from the corresponding symmetric pressure patch. The response from the actual pressure patch that is lost across the symmetry line is identical to the response from the corresponding symmetric pressure patch. This modified height basis function is now symmetric about the centerline and can be safely used while still only using the pressure patches on one side of the hull and computing the wake over the same half of the centerline. 02468100123456xy

Fig. 24 Total Constant Pressure Patch Basis Function

4.2. SINUSOIDAL PRESSURE METHOD SYMMETRY

A similar symmetry condition exists for the sinusoidal pressure method; the symmetry is easier in this case though. The pressure basis functions in this method are defined as existing over the entire width of the hull. This means they can be inherently symmetric without reflecting the response of another pressure patch. From (2.25) the pressure is defined as a series of sine waves with increasing frequency. Essentially the pressure is approximated with a simplified Fourier series. This pressure is symmetric if

75

only odd terms are used for n. To verify this symmetry condition a sinusoidal pressure patch is simulated without the symmetry assumption. As can be seen in Fig. 25, the resulting surface disturbance is indeed symmetric. It is valid to only compute half of the wake. xy0246810-505

Fig. 25 Sinusoidal Pressure Patch Wakes

To strictly satisfy the derivation all Fourier terms are used to define the pressure, but because of symmetry if the hull (and thus pressure disturbance) is symmetric the coefficients corresponding to even n are zero. Computation time can be saved, however, by simply discarding these terms and only computing the response from the odd terms of n. To take this into account the index for the pressure terms is modified.

*2nn 1 =− (4.1)

where . This modified index is always odd. It is used in place of the old index in each of the sinusoidal pressure method basis functions. Now the pressure basis function is guaranteed to be symmetric and the resulting wake will be as well. This means that once again only half the wake has to be simulated. Also, only half the basis *1,2,3,n=L

76

functions need to be determined because each of the asymmetric pressure basis functions is ignored.

4.3. SINUSOIDAL PRESSURE METHOD OSCILLATION

The fact that the frequency is increasing in the sinusoidal pressure method brings up an interesting facet of this method. To accurately capture the characteristics of a sine wave a minimum of two points per period of the sine wave must be known. This is the minimum frequency (Nyquist frequency) at which a signal must be sampled in order to reconstruct it using a Fourier series [Franklin]. In practice, it is desirable to sample at more points per period than the Nyquist frequency suggests because a Fourier series is not generally used to reconstruct the signal; a simpler method (e.g. linear interpolation) is used in reconstruction. These simpler methods perform better with more than two samples per period.

The number of points this method is allowed is limited, however. In order for there to be as many pressure coefficients as there are calculated wave heights under the hull the number of points from the centerline to the beam must equal the number of pressure terms. The number of samples per period is obviously lowest with the highest frequency. The worst case number of samples per period can therefore be calculated by comparing the number of points per half-beam to how many periods there are in a half-beam at the highest frequency. This ratio is

worstcase#21samplesNperiodN= 4

− (4.2) 77

where N is the number of pressure terms. In the limit as more pressure terms are used, the number of samples per period approaches the Nyquist limit.

0510152022.533.54Number of Pressure TermsMinimum Number of Samples /

Period

Fig. 26 Loss of Data as More Fourier Terms are Used

It does not take many pressure terms before there are essentially two samples per period in the highest frequency pressure basis function. This means the higher frequencies are in danger of loosing data due to not sampling often enough. One method to examine the impact of this problem is to oversample the system and compare the results to the results produced by the actual sampling.

The first function to study is the pressure basis function. The lowest and highest frequency pressures are shown along with where they are actually sampled and a linear interpolation (Fig. 27). The maximum number of pressure terms is 10.

78

00.0750.1750.2750.3750.475-1-0.500.51yP* n = 1 – oversampledn = 1 – actualn = 10 – oversampledn = 10 – actual

Fig. 27 Oversampled Pressure Basis Function

It is obvious in this plot that the sampling captures the character of the low frequency pressure, but does not capture the high frequency oscillations. The amplitude of the high frequency oscillations remains constant, but the amplitude of the sampled signal seems to increase moving away from the centerline.

It seems likely that similar behavior will be observed in the wake basis function that results from this pressure. The model does not guarantee that the oscillations carry over to the wake, but physically it makes sense. The plot only shows one slice of the wake, but the behavior is similar throughout the entire wake (Fig. 28). Once again the low frequency component is captured well. A blip at y = 0.175 is not captured, but it is most likely noise and should be ignored. As with the pressure basis function the high frequency component, on the other hand, is not captured. The amplitude actually decreases slightly, but the sampled function has a peak amplitude in the middle of the hull then it decreases again.

79

00.0750.1750.2750.3750.475-0.500.51yz n = 1 – oversampledn = 1 – actualn = 10 – oversampledn = 10 – actual

Fig. 28 Oversampled Wake Basis Function

m = 1, x = -1.1

The data loss can also be examined by looking at the characteristics of the total wake basis function compared to an oversampled basis function (Fig. 29). The oversampled function shows the oscillations are reasonably consistent in both directions; the peaks and valleys create a regular grid. The sampled basis function, however, seems to have a series of nearly plane waves with peaks at a slight angle from the centerline. The two behaviors are quite different. It is also important to note that the data loss issue discussed above is in the transverse direction. A similar behavior is observed in the longitudinal direction.

80

ActualOversampled

Fig. 29 Oversampled Wake Basis Function Surface

m = 1, n = 10

This problem becomes particularly pressing when the hull shape is known and the pressure coefficients are unknown. Because of the data loss the higher frequency basis functions become closer to linearly dependant and the matrix inversion harder to accurately compute. In practice, a high frequency term in a forward pressure strip is erroneously given a large amplitude. This creates large high frequency oscillations further aft along the hull which are fixed by assigning an even larger amplitude to the high frequency terms in these pressure strips. The inversion results in an estimated hull that when oversampled is nothing but an exponentially growing high frequency oscillation. The wake this hull creates is obviously not accurate.

In conclusion, it is possible to take advantage of the symmetry of boat hulls to reduce the number of calculations needed. With the constant pressure method, this requires slightly modifying the basis function. For the sinusoidal pressure method all that is needed is to discard the pressure terms that are asymmetric. These procedures effectively reduce the computation time for each method. The fact that the sinusoidal pressure method is based on a Fourier series also raises issues related to the number of points used

81

to capture the system behavior. In general, these issues along with all the inadequacies discussed in chapter 3 make the sinusoidal pressure method unacceptable for computing the shape of a boat wake. Throughout the remainder of the work the constant pressure method is used exclusively.

4.4. ISSUES WITH INVERSION

As was mentioned in chapter 2, the wave height over a grid of points is determined by summing a series of basis functions multiplied by pressure coefficients at each point (2.170). In other words, the wave heights are a linear combination of basis functions:

gρζ= AE. (4.3)

The ideal method to determine a boat’s wake is to define the surface of the hull and compute the basis functions over the wetted area. The number of points defining the hull should be the same as the number of unknown pressure coefficients (hullE is square). The pressure coefficients are then determined using

()1hullhullAgρ ζ −=E. (4.4)

A new set of basis functions is then computed over the wake area and the pressure coefficients determined with the hull are used to solve for the height of the wake:

.wakewakeAgζρ=E (4.5)

The problem with this method is that it involves inverting the basis function matrix. Because of the uncertainties that go into creating this matrix the inversion is not a stable numerical operation.

82

While the basis function matrix is nonsingular (theoretically invertible) it is very poorly conditioned and thus sensitive to noise. The condition number of a matrix is the ratio of the largest and smallest singular values of the matrix. For the ideal case (an orthogonal matrix) the condition number is one. As the condition number grows the matrix inversion becomes less stable and the accuracy of the inverted matrix decreases. Chapter 3 presents two cases to validate the solution methods, a sinusoidal pressure and a simplified hull pressure. In these cases the condition numbers of the hull basis functions are 7.5·109 and 1.2·1011 respectively. In these tests the pressure (as opposed to the hull shape) was defined; however it would not be safe to invert the hull basis function matrices.

To demonstrate the instability the basis functions from the simplified hull pressure will be used with two hull shapes and the resulting pressures compared. The first hull shape is the one computed using constant pressure patches (Fig. 30).

Fig. 30 Original Test Surface

When the basis function matrix is inverted the pressure (Fig. 31) is identical (to machine precision) to the actual pressure used to generate the surface (Fig. 18).

83

Fig. 31 Pressure Resulting From Original Test Surface

The surface is then rotated by 0.00001° about the y – axis (adjusting the trim by a miniscule amount). This small change in the hull should produce virtually no change in pressure; however, the pressure solved for (Fig. 32) is unrealistic. The pressure rise at the front of the surface displacement is still visible. However, further aft the pressure oscillates uncontrollably. As the surface is rotated more (or other perturbations are applied) the amplitude of the oscillation becomes much larger. 84

Fig. 32 Pressure Resulting From Rotated Test Surface

To better understand the reason for these oscillations the basis function matrix must be examined. The hull shape equation (4.3) can be written in block form as a series of smaller matrices due to its formation from a multidimensional tensor equation, i.e.

1,11,1,1,11,2,1,11,,1,11,12,11,1,2,11,2,2,11,,2,11,2,11,1,,11,2,,11,,,11,hullJNJNJNMJJNJNJNMJNIJNIJNIJNMIJNMgAEEEAEEEAgEEEAρζζζρζ→→→→→→→→→→→→→→→→→→→→→→→→=⎡⎤⎡⎤⎡⎢⎥⎢⎥⎢⎢⎥⎢⎥⎢=⎢⎥⎢⎥⎢⎢⎥⎢⎥⎢⎣⎦⎣⎦⎣ELLMMMOML⎥⎥

N ⎤⎥⎥

⎦

M

(4.6)

where ,ijζ, , and are as in (2.170)and “→” represents the whole range of values (i.e. a vector/matrix). The basis matrix can be simplified by noting that the contribution of a pressure patch only affects the flow aft of the patch: ,,,nmijE,nmA

,,,0 ; 1nmijEi m =≤−. (4.7) 85

This leads to:

86

⎤⎥⎥⎥⎥⎦

. (4.8) 1,1,1,11,11,11,1,2,11,2,2,12,11,2,11,1,1,,11,2,,11,,,100000NJJNNJNJJNIJNMNIJNIJNMIJEAEEAgAEEEζζρζ→→→→→→→→→→→→→→→→→→⎡⎤⎡⎤⎡⎢⎥⎢⎥⎢⎢⎥⎢⎥⎢⎢⎥=⎢⎥⎢⎢⎥⎢⎥⎢⎢⎥⎣⎦⎣⎢⎥⎣⎦LLMMOMMML

The first row of hull heights (1,1Jζ→) can be determined independent of any pressures except for the first row (). Conversely, the first row pressure vector can be found by inverting the first row basis function matrix () to match the first row hull height vector. This pressure creates a disturbance down stream. The second row of pressures must cancel this disturbance as well as match the second row height to the hull. Due to noise in the basis function calculations, the down stream height disturbances contain small oscillations. To cancel the oscillations, the pressure in the next inversion contains the same oscillation with a slightly larger amplitude. As the solution propagates aft the oscillation in the pressure continues to grow causing the instability. 1,NA→ 1

J 1,1,1,1NE→→

0.10.20.30.4-15-10-50510yPx = 1.5x = 1.7x = 1.9

Fig. 33 Increasing Pressure Oscillation

This instability and sensitivity to small hull shape changes require a change in solution procedure. Instead of defining the hull, determining the pressure that produces that hull, and finally determining the wake the new procedure is to define a pressure shape with unknown parameters; vary those parameters over a reasonable range; and then determine both the hull and wake shape the pressure produces. One potential problem with this method is that the hull shapes produced may not be realistic. The hull shape needs to be monitored to guarantee that it is practical. This may also be a benefit, however; since the hull is not being varied directly, some shapes may be discovered that are possible to manufacture but would not otherwise be investigated.

87

88

5. HULL OPTIMIZATION

The goal of this research is to study methodology for designing hulls for different types of professional and recreational tow boats, specifically for towing a skier or a wakeboarder. Algorithms for predicting the wake (and hull shape) given a pressure disturbance have been developed. The next step is to describe what is “good” for a wake and what parameters are defined in each event. The wake shape equations depend on the speed the boat is traveling; the speed is in turn dependent on what and who is being towed. Additionally, it does not make sense to examine the shape of the entire wake; the rider will only be in a portion of the wake defined by the tow rope length which depends on the event and the rider. Different events require different parameters and different riders have different preferences.

5.1. EVENT DEFINITION

Skiing has existed for the longest time and its events are the most formalized. In particular two governing bodies oversee most of the professional skiing tournaments, the International Water Ski Federation (IWSF) and the American Water Ski Association (AWSA). These two bodies have very similar rules concerning boat speed and rope length for events. Two events will be considered for skiing: jumping and slaloming.

Jumping is defined as the skier using a ramp to jump and get the maximum distance before landing. The skier can choose a speed between 45 and 57 kph (28 and 35 mph) and the tow rope length is 23 m (75 ft) [AWSA].

Fig. 34 Ski Jumping Example and Course Layout

Slaloming is defined as the skier going around a series of buoys to either side of the boat’s path. The speed is set at 58 kph (36 mph); the rope length varies between 11.25 and 18.25 m (37 and 60 ft) [AWSA]. The rope gets shorter with each run making the event more difficult.

89

Fig. 35 Slaloming Example and Course Layout

The conditions for recreational (amateur) waterskiing vary much more based on the skier’s preference. Typically the speed is less than for professionals and the tow rope length is towards the longer end of the professional range. One factor ties each of these activities together: the activity is not based on the wake. The boat must be towing the skier, but fewer disturbances from the wake are better.

Wakeboarding is a newer sport; two bodies also govern its professional events, the World Wakeboard Association (WWA) and the World Wakeboard Council (WWC). Tournaments sponsored by these two bodies have different events that mainly consist of various freeform jumps, flips, etc. Both agree, however, that boat speed and rope length are at the discretion of the rider [WWC]; in general the riders prefer boat speeds much slower than skiing, 12 to 20 mph (19 to 32 kph), and tow rope length in the middle of the slalom lengths [Favret, Solomon]. This is also true for recreational wakeboarding. One factor that is consistent for most wakeboarding activities is that the wake itself provides the launching platform for many of the tricks. Therefore, the wake should be larger and provide a ramp for the wakeboard.

90

Fig. 36 Wakeboarding Examples

In summary, besides the actual hull shape there are two parameters that are varied for the different activities: the boat speed and the tow rope length, or distance behind the boat the wake will be examined. Two basic types of wake are also desired depending on the activity. For skiing a flat wake as small as possible is desired. For wakeboarding a larger wake with ramp characteristics is ideal.

To characterize these wake shapes two parameters will be examined: the maximum height of the wake and the slope of the main disturbance. These parameters are calculated along the arc the skier makes behind the boat and are defined as

()()1tan.maxminmaxminmaxminheightzzzzslopeRzRzθθ−=−−⎡⎤=⎢⎥−⎣⎦ (5.1) 91

z

θ

zmin

zmax

()zmaxθ (z )

θ min

R

θ

y

x

Fig. 37 Wake Parameters

For each speed and tow length under consideration the pressure distribution that creates the hull will be varied to attempt to either minimize or maximize (depending on the activity) these parameters. The hull shape that created the wake will then be examined for reasonability and shape.

5.2. TYPICAL HULL SHAPE

Because the hull shape is being determined from the surface deformation made by the pressure instead of being defined directly, the shape must be examined to guarantee that it is reasonable. The gross shape of most planing hulls has not changed drastically in

92

many years, although minor improvements have been made. The shape is characterized by a “V”-bottom that is sharper near the bow with decreasing deadrise angle towards the stern. This basic shape is shown below [Comstock]. LCLB0 ½ 1½ 1 2 3 4 5 6 7 8 9 10

Fig. 38 Basic Planing Hull Shape

Two critical parameters characterize this shape. The chine is the sharp break in the transverse slope of the hull (the edge of the “V”-bottom). The slope of the “V” is known as the deadrise; as was previously mentioned it varies along the length of the hull.

A third parameter needed to characterize the hull’s shape in the water is the trim angle. This angle is how parallel the keel is to the free surface. This is not a parameter of the hull itself; the trim angle is maintained to balance the pitch moments acting on the hull.

Modern planing hulls are more complicated than this example. For instance most have strakes, longitudinal ridges running along the bottom of the hull; however, the chine, deadrise, and trim angle characterize the overall free surface deformation created by the hull well enough for the purposes of this research.

93

94

5.3. TEST CASES

As mentioned in section 5.1, the conditions that represent each event include speed and tow rope length. Six cases are studied to find optimal hull shapes. To represent ski jumping the rope length is defined as 23 m and the speed can vary between 45 and 57 kph; these limits are tested. In slaloming the speed is set at 58 kph and the rope length can vary between 11.25 and 18.25 m; 11.25 and 14.25 m are evaluated. Wakeboarding is not regulated but two typical cases, 14.75 m rope length at 19 kph and 32 kph, are tested. These parameters are summarized in Table 1.

Table 1 Test Parameters

Speed (kph)

Tow Rope Length (m)

Jump 1

57

23

Jump 2

45

23

Slalom 1

58

11.25

Slalom 2

58

14.25

Wakeboard 1

19

14.75

Wakeboard 2

32

14.75

Three other parameters must be set for each test: the displacement, length, and beam of the boat. When planing the volume of water displaced does not balance the weight of the boat due to the lift created by planing. Instead of measuring displacement the lift generated by the pressure is adjusted to match the weight of the boat. The length and beam are more difficult. As more lift is created the boat rises out of the water. Therefore, the wetted area will not necessarily stay the same as parameters change. In

order to minimize variations and concentrate on how hull shape (as opposed to size) affects the wake for this study the length and beam remain constant for each test.

From section 5.2, planing hulls normally have a chine where the shape drastically changes; above the chine the hull is nearly vertical. When planing the waterline is very near the chine (see Fig. 39). Therefore, the wetted beam is assumed to be 90% of the overall beam. Similarly, approximately two thirds of the length of the hull is in the water when planing (see Fig. 39). For the tests a typical tow boat with parameters in Table 2 is used.

Fig. 39 Examples of Planing Boats

95

Table 2 Typical Tow Boat Parameters

Overall Value

Wetted Value

Mass (kg)

1445

1445

Length (m)

6.45

4.3

Beam (m)

2.31

2

The pressure distribution is varied for each test. The pressure distribution used in section 3.2 represents a simplified pressure produced by a planing hull. In particular it has characteristics that are present in planing craft such as a singularity at the leading edge and zero value at the other edges [Scullen]. These characteristics can be maintained and the pressure allowed to vary by modifying the pressure equation to

11cdLyPxB⎡⎤⎡⎤⎛⎞⎛⎞=−−⎢⎥⎢⎥⎜⎟⎜⎟⎝⎠⎝⎠⎣⎦⎣⎦. (5.2)

In this equation c effectively controls how long the peak of the pressure is and d controls how wide the peak is. For the skiing tests (jumping and slaloming) the goal is to minimize the wake height and slope; for the wakeboarding tests the goal is opposite, to maximize the wake height and slope.

In each test using the conditions listed above (Table 1 and Table 2) the pressure is varied according to equation (5.2) and then modified so that lift matches weight. The wake generated by this pressure is computed and examined in an arc defined by the tow rope. The maximum wake height and wake slope along this arc are determined using equation (5.1). These two parameters are then compared over the range of pressure parameters so that for each test the slope and height can either be maximized or minimized as appropriate.

96

97

5.4. JUMP 1 TEST RESULTS

The first test is the higher speed ski jump test. The surfaces shown below represent the maximum wave height and slope as a function of the two pressure parameters (c and d). The goal of this test is to minimize both the maximum height and the slope. The maximum height is minimized when either c or d is at the minimum value tested; the absolute minimum occurs when both c and d are at the minimum tested value. The minimum value of c tested is 0.01. Values below this do not appreciably modify the pressure from this value and so are not tested. Similarly for d testing any lower does not change the results.

The minimum slope occurs where d is approximately 4 and c is greater than 5. It stays near this minimum along the line c = 3.5. To compromise between these two minimizations values of c = 2.3 and d = 2 are chosen for the optimal wake.

c d246810510152025300.40.450.5 c d246810510152025302829303132

Fig. 40 Jump 1 Minimization

5.4.1. Resulting Pressure and Hull Shape

These pressure parameters result in the pressure shown below (Fig. 41). This pressure produces the wake also shown below. In general this hull shape looks reasonable. The trim angle is 5.7 deg; at the stern of the boat the deadrise is 6 deg. Unlike the typical hull presented in section 5.2 the deadrise decreases toward the front of the boat. This is due to two reasons. First, the wetted area is defined as a rectangle instead of the triangle that a more typical hull would have. Without knowing the hull shape beforehand it is not possible to know the shape of the wetted area triangle so a

98

99

rectangle must be assumed. Second, near the bow the pressure function used (5.2) has value for the width of the boat instead of only near the centerline (again a rectangular wetted area instead of a triangle). A different pressure function might produce a more typical hull shape. One interesting feature of the hull shape is the positive height on the outside rear corner. At the transom approximately 0.8 m from the centerline the hull begins to rise quickly. This is characteristic of the chine line. The chine widens as it approaches the transom; this is not typical, but has the same explanation as the atypical deadrise. Also, outside the chine the height becomes positive. It is common for the water to spray up on the sides of the boat (see Fig. 39). The positive height is physically reasonable in this light.

x (m) y (m)12340.20.40.60.8123456x 104 x (m) y (m)12340.20.40.60.81-0.4-0.200.2

Fig. 41 Jump 1 Pressure and Hull Shape

5.4.2. Resulting Wake

The wake this pressure produces is shown below in Fig. 42. The black line represents the arc a skier would follow. The height the skier encounters is also plotted vs. the distance along the arc. The square indicates the highest point on the arc and the circle the lowest. A typical wake behind a boat is reasonably flat directly behind the boat with a crest spreading in a “V” outside this flat area. Outside the crest is a valley which levels off to the undisturbed free surface height. The shape shown below has several of these characteristics. In particular the size and location of the crest and valley are the

100

101

parameters used to characterize the wake. The area inside the wake crest (inside the “V” behind the boat) is not accurate. This area has longitudinal oscillations that are not present in physical wakes. These oscillations narrow further behind the boat until the surface inside the wake becomes essentially noise (around x = 25 m in this test). These oscillations can also be seen in the centerline validation tests (section 3.3). However, in this area of the wake the extra disturbance from the boat’s propeller affects the fluid flow; the disturbance from the propeller is not included in this simulation. Therefore, the area inside the “V” behind the boat can safely be ignored when analyzing the wake shape.

Ignoring the area directly behind the boat, the height profile along the skier’s arc appears reasonable. The valley outside the wake is not typically that much more extreme than the surface height inside the wake crest, but this is most likely due to a combination of the atypical hull and the no propeller assumption. The height of the wake is 39.8 cm and the slope is 28.4 deg.

x (m) y (m)5101520251234-0.4-0.200.2 12345-0.3-0.2-0.100.10.2Distance Along Arc (m) z (m)

Fig. 42 Jump 1 Wake

5.5. JUMP 2 TEST RESULTS

The next test case is the slower ski jump. As before, the surfaces to minimize (wake height and slope) are shown below. The height shows a strong trend: as c is decreased the height decreases. The slope is not as clear. There is a rough minimum along the line c = 2.1 and the line d = 2 is also a minimum for c greater than 4. The slope is also near minimum at c = 0.01 and d = 14.6; this is the condition chosen as the optimum compromise. 102

c d246810510152025300.560.580.60.620.640.66 c d246810510152025302525.52626.527

Fig. 43 Jump 2 Minimization

5.5.1. Resulting Pressure and Hull Shape

These parameters result in a pressure that is nearly linear while still maintaining the required boundary conditions. The hull shape is also interesting. The trim is normal at 5.2 deg, and as before the hull gets wider forward. The deadrise, however, is not apparent. In fact, the rear outside edge is lower than the centerline. This shape is known as a catamaran and while not typically used in towboats is common in other power boat designs.

103

x (m) y (m)12340.20.40.60.81000200030004000500060007000 x (m) y (m)12340.20.40.60.81-0.4-0.200.2

Fig. 44 Jump 2 Pressure and Hull Shape

Although this hull shape is not common for towboats, some manufacturers are beginning to produce boats with similar features. The boat shown below is the MasterCraft X-Star. It is not a true catamaran because the “V” hull is still present; however, outriggers are present similar to a tri-hull. The X-Star is a wakeboarding tow boat; it is not meant to imply that the hull shape presented above is necessarily the best ski tow boat. It does, however, show that similar designs are feasible and being used in industry.

104

Fig. 45 MasterCraft X-Star

5.5.2. Resulting Wake

As before, the wake looks reasonable except for the area immediately behind the boat. The maximum wake height is 56.2 cm and the slope is 25.9 deg. In general, this wake is not as preferable as the previous test case. As boat speed increases typically the wake becomes smaller and flatter because less of the hull is disturbing the free surface to achieve the same lift. The comparison of the two speeds qualitatively matches the physical expectations. The magnitude of the height for both this and the previous test seems high. More evaluation of this technique is warranted to guarantee its accuracy.

105

x (m) y (m)5101520251234-0.500.5 12345-0.2-0.100.10.2Distance Along Arc (m) z (m)

Fig. 46 Jump 2 Wake

5.6. SLALOM 1 TEST RESULTS

The next two tests are the slalom runs. The first slalom test is with the shorter rope. The trends in both the height and slope are strong. The height decreases as either c or d decrease; the absolute minimum occurs at c = 0.01 and d = 2. The slope has similar trends except the minimum in c occurs at 2; for c < 2 the slope begins to rise again. The gradient in the d – direction is stronger, however (d has more affect on slope than c). Because c = 0.01 and d = 2 minimizes the height and is near the minimum slope this point is chosen as the optimal compromise.

106

c d246810510152025300.40.450.50.550.6 c d246810510152025302829303132

Fig. 47 Slalom 1 Minimization

5.6.1. Resulting Pressure and Hull Shape

The pressure this condition produces is fairly linear longitudinally and curved transversely. The hull, however, does not exhibit the catamaran shape of the previous linear hull; it is a more traditional “V” shape. The trim angle is 4.9 deg and the deadrise at the transom is 2.6 deg.

107

x (m) y (m)12340.20.40.60.82000400060008000 x (m) y (m)12340.20.40.60.81-0.3-0.2-0.100.10.2

Fig. 48 Slalom 1 Pressure and Hull Shape

5.6.2. Resulting Wake

The wake this pressure produces is slightly different from the previous two test cases. The lowest height of the wake occurs inside the wake crest as opposed to the valley outside the crest. In fact, no valley outside the wake crest is present. This difference is due to the fact that the skier is so close to the boat; further from the boat the wake crest for this peak dies and the wake crest the previous tests used is formed. This close to the boat the skier is ahead of the formation of the main wake. The height the skier sees is 38.4 cm and the slope is 28.2 deg.

108

x (m) y (m)5101520251234-0.4-0.200.2 12345-0.4-0.3-0.2-0.100.1Distance Along Arc (m) z (m)

Fig. 49 Slalom 1 Wake

5.7. SLALOM 2 TEST RESULTS

The longer rope slalom test results are interesting because for many pressure distributions the wake is not valid; in particular the main wake crest is below the undisturbed free surface which is not realistic. These obviously incorrect wakes are not considered in the optimization. Of the wakes that are reasonable the heights and slopes follow the trends seen previously in the minimization. As c and d decrease the height and slope also decrease; d has a bigger effect on slope than c does. The optimal combination for this test is c = 2.3 and d = 17.4.

109

c d246810510152025300.50.520.540.560.58 c d246810510152025303939.54040.5

Fig. 50 Slalom 2 Minimization

5.7.1. Resulting Pressure and Hull Shape

These pressure coefficients result in a distribution that is wide, but concentrated at the front of the wetted area. The hull has the catamaran shape seen in a previous test. The trim angle is 3.8 deg; the deadrise angle is meaningless in a catamaran hull.

110

x (m) y (m)12340.20.40.60.81234x 104 x (m) y (m)12340.20.40.60.81-0.4-0.200.20.4

Fig. 51 Slalom 2 Pressure and Hull Shape

5.7.2. Resulting Wake

The longer rope in this test means the skier is back in the main wake so the valley outside the wake crest is apparent. The wake height is 48.5 cm and the slope is 39.2 deg. As with the ski jumping tests, the wake heights are excessive compared to physical wakes.

111

x (m) y (m)5101520251234-0.500.5 12345-0.4-0.3-0.2-0.100.1Distance Along Arc (m) z (m)

Fig. 52 Slalom 2 Wake

5.8. WAKEBOARD 1 TEST RESULTS

The last two tests are for wakeboarding. The main difference between these tests and the skiing tests is the lower speeds. Additionally, instead of minimizing the wake height and slope the goal of the tests is to maximize them. The first wakeboard test is the slower of the two. The height and slope show similar trends as the other tests except c has a much larger effect on both the slope and the height than d; in fact unlike before as d decreases the height and slope may increase. The maximum height and slope both occur at c = 10 and d = 30.

112

c d246810510152025301.21.41.61.8 c d246810510152025302426283032

Fig. 53 Wakeboard 1 Maximization

5.8.1. Resulting Pressure and Hull Shape

This pressure is virtually constant over the width of the boat and is concentrated exclusively at the leading edge of the wetted area. The hull shape this produces is problematic. Before the transom of the boat the hull is rising implying the hull is bowl shaped, which is not a practical planing hull design. Considering only the centerline forward of the minimum hull height the trim angle is 35.2 degrees, very steep for a planing hull. As boat speed increases the boat transitions from operating in a displacement mode (lift comes from displaced water) to a planing mode. During this transition the boat trim is very large. The hull in this test with these pressure parameters 113

is going slow enough that it has not made the transition to full planing. The wake shape equations are only valid for full planing. This is the most likely cause for the unrealistic hull shape. x (m) y (m)12340.20.40.60.82468x 104 x (m) y (m)12340.20.40.60.81-1.5-1-0.50

Fig. 54 Wakeboard 1 Pressure and Hull Shape

5.8.2. Resulting Non-Optimal Pressure and Hull Shape

To find a reasonable hull the parameters which minimize the wake height and slope (c = 0.01, d = 2) are examined. The hull this produces still curves upward at the stern, but not much.

114

x (m) y (m)12340.20.40.60.82000400060008000 x (m) y (m)12340.20.40.60.81-1.5-1-0.50

Fig. 55 Wakeboard 1 Non-Optimal Pressure and Hull Shape

5.8.3. Resulting Non-Optimal Wake

The wake the non-optimal pressure coefficients produce is obviously not the best wake, but it does look reasonable. The crest “V” is very wide which is typical of slower boats, particularly when not fully transitioned to planing mode. The lowest point of the wake is inside the crest which is not unreasonable. The wake height is 110.6 cm and the slope is 28.5 deg. As with all the other tests this height seems excessively large.

115

x (m) y (m)5101520251234-1.5-1-0.500.5 1234500.20.40.60.8Distance Along Arc (m) z (m)

Fig. 56 Wakeboard 1 Non-Optimal Wake

5.9. WAKEBOARD 2 TEST RESULTS

The final test is the faster wakeboard test. The height and slope profiles have similar trends as the slower wakeboard test. The maximum of each occurs at c = 10 and d = 30.

116

c d246810510152025300.650.70.750.80.850.90.95 c d2468105101520253031.53232.53333.534

Fig. 57 Wakeboard 2 Maximization

5.9.1. Resulting Pressure and Hull Shape

This pressure is identical to the optimal wake in the previous wakeboard test, but this hull is more reasonable. One interesting feature of the hull is a depression near the forward outside corner of the wetted area. This is similar to a catamaran, but vanishes at the stern of the boat. No tow boats are currently manufactured in this shape, but it is not infeasible. The trim angle is 11 deg and there is no deadrise. 117

x (m) y (m)12340.20.40.60.82468x 104 x (m) y (m)12340.20.40.60.81-1-0.500.5

Fig. 58 Wakeboard 2 Resulting Pressure and Hull Shape

5.9.2. Resulting Wake

This test has a high enough speed that the wake looks more like the skiing wakes. The valley outside the wake crest is the lowest point although the valley inside the crest is not much higher. The wake height is 95.3 cm and the slope is 34.3 deg. It does not seem reasonable for a boat to create a wake nearly one meter high.

118

x (m) y (m)5101520251234-1-0.500.51 12345-0.6-0.4-0.200.20.4Distance Along Arc (m) z (m)

Fig. 59 Wakeboard 2 Wake

5.10. RESULTS SUMMARY

In each test the wake the boat produced appears reasonable. A “V” shaped crest is visible with a valley inside it and outside either a valley or drop to the undisturbed free surface. The wake is wide when the boat is slow and narrows as speed increases matching reality. The area directly behind the boat is not realistic, but this could be due to ignoring propeller effects. The wake shape also shows a strong correlation to the pressure distribution, and the hull shapes while not always typical are usually reasonable.

There is however a problem with the wake height. In every case tested the height is larger than can be expected in a real wake. The shape of the wetted area also needs to be

119

120

examined. In particular, each hull produced is flat at the leading edge of the wetted area instead of showing the characteristic deadrise present in tow boats. Without the ability to directly define the hull shape defining the shape of the wetted area becomes problematic.

121

6. CONCLUSIONS AND SUGGESTIONS FOR FUTURE RESEARCH

6.1. CONCLUSIONS

In conclusion, this work attempts to predict the wake shape behind a given boat moving at a given speed. To do this two methods for determining the free surface deformation from a pressure distribution are presented. Both methods derive from the same velocity potential function; the difference is the form of the pressure distribution and the resulting integrals that must be evaluated numerically. Each method is compared to previously published results to verify its accuracy. The sinusoidal pressure basis function method does not yield results matching the original authors. Furthermore, this method has potential instabilities beyond the numerical inaccuracies; in order to produce a square linear system the pressure oscillates too quickly for the solution grid to capture. The constant pressure basis function method yields results that closely match the published results. The numeric algorithm for the constant pressure basis function also runs much faster than the sinusoidal pressure basis function. Because of the sinusoidal pressure basis function method’s shortfalls, the constant pressure basis function method is used to predict the wake shape behind the boat. Even the constant pressure basis function method is not stable enough to invert and determine a pressure distribution from a boat hull shape. Instead a variety of pressures are examined for wake shape and hull shape.

122

Tests are performed for a range of scenarios encountered in real life. In particular, ski jumping and slaloming, and wakeboarding conditions are simulated. To find the optimal hull shape for each of these activities the wake shapes produced by the various pressure distributions are examined for characteristics important to the activities: wake height and slope. The hulls produced by the tested pressures are feasible for the most part. However, because the hull is not directly defined some characteristics are not typical. Specifically, every hull produced is flat at the leading edge of the wetted area; this is not typical of modern towboats. The flat hull shape is due to the fact that a rectangular wetted area must be assumed. Additionally, several hulls had depressions along the outside edge. This shape is common in power boats, but not in tow boats. The wakes observed in the simulations are also qualitatively reasonable. A definite “V” crest is visible behind the boat. This “V” widens as the boat moves slower matching physical intuition. The wake height and slope also show a strong correlation to pressure distribution. Oscillations appear inside the “V” that do not represent reality, but these are most likely due to ignoring the effects of the propeller. The major discrepancy between the simulated wake and a physical wake is the difference in height. The simulated wakes are much larger than is typical of tow boat wakes. This discrepancy warrants further research.

One final important note is that the hull optimization performed in this research only considers wake parameters. Boat hulls must be designed with many other considerations such as handling, drag (engine size), and ride comfort in mind. For example, most of the hulls studied in this research have nearly flat bottoms. Flat bottom boats are common for small one or two person, low power boats for use in calm waters. Larger boats (such as

123

tow boats) never have flat bottoms because of considerations completely outside wake shape. Flat bottom boats are not as maneuverable as “V” hulls; they slide over the water instead of cutting through it. This sliding behavior also makes their ride in rough water very violent; the boat goes up over each wave instead of cutting through them. Any designer who wished to make use of the procedures developed in this research would need to combine the wake considerations with other boat performance considerations. Until a method to invert the surface deformation procedure (define a hull shape instead of a pressure distribution) is developed it will be very difficult to use this method in conjunction with other hull considerations.

6.2. SUGGESTIONS FOR FUTURE RESEARCH

There are several areas that could be investigated to improve the method developed in this research. These areas generally fall into two categories: numerical algorithms used and range of tests to optimize hull shape. The most complex algorithms in this research involve the numerical integration of functions. These functions have features making them difficult to integrate, singularities and rapid oscillations. Transformations are used to remove the singularity and slow the oscillations, but these change the integration range to negative and positive infinity introducing inaccuracies when the integration is truncated. Better analytical transformations may be possible to improve the numerical integration accuracy. The particular pressure basis functions used can also improve the numerical algorithms. The sinusoidal pressure basis function method has the problem of data loss from not capturing higher frequency oscillations. The method does satisfy the Nyquist criterion for number of points per period necessary to uniquely define a sine wave; however, this absolute minimum is not enough to accurately capture the sine wave.

124

This apparent discrepancy is due to the fact that the Nyquist criterion is for fitting the data with a Fourier series instead of the piecewise linear interpolation used in the method developed for this research. If the wake surface were described by a Fourier series instead of being determined by inverting the basis function matrix it is possible the higher frequency oscillations could be captured thereby stabilizing the sinusoidal pressure basis function method. Another improvement that could be made with the pressure basis functions is in the constant pressure basis function method. This method approximates the pressure distribution by piecewise constant rectangular patches. This creates discontinuities in the applied pressure disturbance. The free surface disturbance may be smoother if these pressure discontinuities were not present. Specifically if the pressure distribution could be approximated by either two dimensional piecewise linear or cubic spline rectangles a smoother wake may result. For this method to succeed the velocity potential function would need to be derived for this new pressure type. All of the numerical improvements have the same goal in mind: to allow the free surface deformation problem to be inverted. For this research to move from an analysis tool to a design tool, a way to define the hull shape, and thereby determine the wake shape, must be developed. In particular this means lowering the condition number of the surface basis function matrix. The condition number may be lowered either by more accurate numerical integration or by using better pressure basis functions.

The second category for recommended future research is in the range of tests to optimize the hull. If a method to directly define the hull shape is developed many more tests become available. One of the more interesting scenarios to investigate is the mixed design tow boat, one that is designed to tow both skiers and wakeboarders. These

125

activities have different wake requirements, but are also performed at different speeds. Could one hull be designed to produce an ideal wakeboard wake at low speeds and an ideal ski wake at high speeds? If the pressure must be defined instead of the hull shape then a wider range of pressure distributions should be examined. Specifically the knowledge that a “V” hull produces a triangular wetted area should be taken into account. Regardless of whether the hull shape or pressure distribution is defined the results need to be compared to physical experiments. This work highlighted results that do not seem physically reasonable. These results need to be compared to a real wake.

126

REFERENCES

1. Abramowitz M. and Stegun I. A., Handbook of Mathematical Functions, New York, Dover, 1972.

2. AWSA, AWSA Official Tournament Rules, 22 July 2005, American Water Ski Association, 10 Oct. 2005, <http://www.usawaterski.org/pages/divisions/3event/ 2005AWSARulebook.pdf>.

3. Comstock J. P. ed., Principles of Naval Architecture, New York NY, The Society of Naval Architects and Marine Engineers, 1967.

4. Cheng X. and Wellicome J. F., “Study of Planing Hydrodynamics Using Strips of Transversely Variable Pressure”, Journal of Ship Research, vol. 38, 1994, pp 30-41.

5. Currie I. G., Fundamental Mechanics of Fluids, 2nd Ed, New York NY, McGraw-Hill., 1993

6. Doctors L. J., “Representation of Three-Dimensional Planing Surfaces by Finite Elements”, Proceedings of the 1st Conference on Numerical Ship Hydrodynamics, 1975, pp 517-537.

7. Favret B. and Benzel D., Complete Guide to Water Skiing, Champaign IL, Human Kinetics, 1997.

127

8. Franklin G. F., Powell J. D., and Workman M., Digital Control of Dynamic Systems, 3rd Ed, Menlo Park CA, Addison Wesley Logman Inc., 1998.

9. IWSF, International Water Ski Federation 2005 Tournament Water Ski Rules, 6 April 2005, International Water Ski Federation, 10 Oct. 2005, <http://www.iwsf.com/ rules2004/ rulebook2004/rules05v1.0.htm>.

10. Lai C. and Troesch A. W., “A Vortex Lattice Method for High Speed Planing”, International Journal for Numerical Methods in Fluids, Vol. 22, 1996, pp 495-513.

11. Maruo H., “Two Dimensional Theory of the Hydroplane”, Proceedings of the 1st Japan National Congress on Applied Mechanics, 1951, pp 409-415.

12. Maruo, H., “High- and Low- Aspect Ratio Approximation of Planing Surfaces”, Schiffstechnik, vol. 72, 1967, pp 57-64.

13. Savitsky, D., “Hydrodynamic Design of Planing Hulls”, SNAME Transactions, 1964, pp 71-95.

14. Scullen D. C. and Tuck E. O., “Free-Surface Elevation due to Moving Pressure Distributions in Three Dimensions”, Journal of Ship Research, 2002.

15. Shen Y. T. and Ogilvie T. F., “Nonlinear Hydrodynamic Theory for Finite-Span Planing Surfaces”, Journal of Ship Research, vol. 16, 1972, pp 3-20.

16. Solomon M., Waterskiing Getting off the Ground, Boston, Aquatics Unlimited, 1997.

17. Sottorf, W., “Experiments with Planing Surfaces”, National Advisory Committee for Aeronautics Technical Memorandum, No. 739, 1934.

128

18. Tong T., A Finite Element Approach to the Planing Problem of High Speed Craft, Diss. University of Southampton, 1990.

19. Wang D. P. and Rispin P., “Three – Dimensional Planing at High Froude Number”, Journal of Ship Research, vol. 15, 1971, pp 221-230.

20. Wehausen J. V. and Laitone E. V., “Surface Waves”, Encyclopedia of Physics, Vol. 9 ed. Flügge S., Berlin, Springer-Verlag, 1960, pp 446-814.

21. Wellicome J. F. and Jahangeer J. M., “The prediction of Pressure Loads on Planing Hulls in Calm Water”, Royal Institute of Naval Architects Transactions, 1978, pp 53-70.

22. WWA, WWA Rule Book, World Wakeboard Association, 10 Oct. 2005, <http://www.thewwa.com/rules.shtml>.

23. WWC, 2005 World Wakeboard Rules, 9 August 2005, World Wakeboard Council, 10 Oct. 2005, <http://www.worldwakeboardcouncil.com/2005worldrules.pdf>.