Journal of Ship Research, Vol. 38, No. 2, June 1994, pp. 9 7 - 1 0 3
II ime°D@maini Caleylatsoris ©f B r i f t Forees ©n FteatBocj
Tw©°yim(grjisö©fiiall O b p e t Bo (Gyr&Brrü airid Wawes
H. J . Prins^ and A. J . Hermans^
in this paper the influence of a current on the drift forces acting on a floating body is studied. A theory has been developed which is correct up to second-order in the current velocity. A numerical algohthm has been found which solves the time-dependent equations and is stable for all current velocities. To test this algohthm, a two-dimensional case has been studied, i.e., the drift forces on a honzontal cylinder of infinite length. A computer program has been developed and very encouraging results have been found for the added-mass and damping coefficients and the drift forces. It is shown that especially for higher velocities the second-order terms in the velocity become very important.
Introduction
I N THE LAST decade numerous attempts have been made to solve the unsteady ship-motion problem. This problem is very important i n predicting the behavior of a ship i n high seas, which includes the interaction between waves and velocity of the ship.
Most of the studies were done using the frequency domain instead of the physical time domain. Unfortunately this re-stricts the theories to harmonic functions. This cannot be realistic of course, for sea waves are not harmonic, especially not i n the neighborhood of a ship. Therefore one has to study the problem i n the time domain. A disadvantage is, of course, the computer time needed to solve the time-dependent prob-lem. Fortunately these computer times become smaller every year with the introduction of new and faster computers. So in the future this disadvantage w i l l not be of importance anymore. Therefore we made a first step i n the direction of solving the time-dependent ship-wave problem. I n order to be able to predict swell i n three dimensions, we first have to look at a two-dimensional (2D) problem. This problem gives us many interesting mathematical and numerical problems which are much easier to solve than i n the 3D problem. So the 2D problem is merely tackled i n order to get some i n -sights into the problems which arise when solving the 3D problem.
I n this paper a horizontal cylinder of infinite length is con-sidered. This gives us a 2D problem which is described i n the first section. Then a numerical algorithm is given which solves the problem i n the time domain correct up to second-order i n the velocity. To test our algorithm we used har-monic functions for the incoming wave field. This enables us to compare our results with the literature. The results found for the added-mass and damping coefficients and the hori-zontal d r i f t force, especially for the lower velocities, agree well with the results found by Zhao (1988), who used the frequency domain instead of the time domain. For higher velocities we found better results due to the fact that we used a second-order free-surface condition instead of the first-order condition used by Zhao. These results are very encouraging to tackle the three-dimensional problem.
'Department of Technical Mathematics and Informatics, Delft Uni-versity of Technology, Delft, The Netherlands.
Manuscript received at S N A M E headquarters December 21, 1992; fi-nal revised manuscript received March 7, 1994.
Mathematical model
We consider a horizontal circular cylinder w i t h radius R (diameter d) and of infinite length floating i n water of infinite depth. The cylinder has draft D. Because of the i n f i -nite length the problem is reduced to a 2D problem. A uni-form current w i t h velocity U and regular incoming waves are traveling i n the positive .x-direction: see Fig. 1. The cyl-inder is free to oscillate i n the x- and 2-directions. The model described below can easily be extended to the case that the cylinder is free to roll as well. The coordinate system is cho-sen i n such a way that the undisturbed free surface coincides w i t h the line z = 0 and the center of the circle is on the line
X = 0. We assume the fluid to be incompressible and we
ne-glect the effect of viscosity. I f we assume the flow to be ir-rotational, we can introduce a velocity potential <I>, which has to satisfy the Laplace equation
V ' * = 0 (1) A t the free surface we combine the well-known dynamic and
kinematic boundary conditions to get
a^<P 1 — + 2V<Ï) • V — - V $ • V(V4» • V<I>) dt^ dt 2 a * - f ^ — = 0 atz=^ (2) dz
w i t h the wave elevation I, given by 1 f l a * '
5 = — - |V<I>f + — \ at free surface
On the ship hull we have the condition that the normal ve-locity of the fluid should match the normal veve-locity of the ship:
— = y„ at instantaneous position of ship hull (3)
èn
To ensure the uniqueness of the solution we have to impose a radiation condition which w i l l be specified later i n this paper.
The velocity potential which satisfies these equations con-sists of a steady and an unsteady part, so we split i t up as
+ + + +
Fig. 1 Geometry
* ( x , t) = <^(x) + 4>(x, t)
Steady potential
In order to calculate the unsteady potential we have to know the steady potential. Unfortunately i t is not well understood how to calculate the wave-making part in the 2D case, for one has to give special conditions at the intersection of the free surface and the body boundary to be able to de-scribe the waves created by the body uniquely. I n the 3D case, however, computation of the steady potential can be carried out by means of a Dawson approach; see for instance Kaven (1992). Therefore we approximate the steady poten-tial by assuming that the steady waves generated by the body are small compared with the incoming waves, which allows us to neglect the wave-making part of the steady po-tential.
The steady potential can then be approximated by ^ ( x ) = Ux
This means that the steady flow is undisturbed by the pres-ence of the cylinder. Of course this is unrealistic, but i t is a first approximation. A better approximation is the double-body potential, which arises from the case when the free sur-face is approximated by a rigid wall, i.e., (d^/dz)/(x, 0) = 0, and when there is no normal velocity at the cylinder. So this potential represents the flow past a lens and is given as a complex function by ^(x) = 2Ua + 1 , _ P c ^ - l (5) with and X + iz
and p and a as i n Fig. 2; see also Moretti (1964). I t can be seen easily that
In the following approximation (5) w i l l be used. Unsteady potential
The unsteady potential consists of all time-dependent parts of the total potential and is therefore composed of the i n -coming wave, the diffracted wave and the potential due to the movement of the ship. So we decompose the unsteady potential as
(\>(X, t) = (}),„,(X, t) + (|)rf,y(x, t) + 4 ) „ ( x , t) (6)
I f we linearize the free-surface condition (2) for fj)(x, t) by assuming ^ to be small compared w i t h cj), which i n t u r n is small, we get dx dxdt 5cj) a (j) 3c|) 3 ~ -F ax' 924, 1 d% + g dx" dz g d z d f j 0 at z = 0 (7) This condition has to be satisfied by or, i f we restrain the cylinder from moving, by 4>inc +
^dif-The well-known Neumann-Kelvin free-surface condition can be obtained by making use of (4) instead of (5) as an approximation for the steady potential:
Nomenclature,
A l l bold symbols represent vectors or matrices.
A = added-mass matrix R = restoring matrix Ir = linearized relative wave height
B = added-damping matrix R = radius of cylinder ^ = second coordinate of Green's
Ci = phase velocity S = computational domain function
D = matrix as = boundary of computational do- p = density of fluid
d = diameter of cylinder main T = dimensionless parameter
D = draft of cylinder t = time (j), <}) = unsteady potential
E = matrix M = size of a time step 4>„ = normal derivative of unsteady
Fn = Froude number U = undisturbed horizontal velocity potential
F l = first-order force of fluid <l>i;ic = potential due to incoming wave
F2 = second-order force X = displacement of cylinder 4i,„ = potential due to movement of
g = gi'avitational acceleration X, z = coordinate system fixed to the cylinder
G = Green's function ship ^dir = potential due to diffracted wave
H = hull surface X = vector i n fixed coordinate system (jj = steady potential
k = wave number
M = mass matrix Greek symbols
(|)d6 = double body potential <E> = total potential (scalar!) n = normal vector, pointing out of 1 = total wave height CO = frequency of encounter
- a
Fig. 2 Lens used for double-body potential
+
2U-dxdt + U a ^ dx'' dz
0 at z = 0 (8) The results obtained by this free-surface condition are clearly inferior to the results obtained by using (5) i n (7).
Furthermore, the linearized version of the body boundary condition (3) becomes
dn
a x
dt n - X • (n • ¥)V<Ï) (9)
with X the displacement vector of the ship, and n the normal pointing out of the fluid domain; see also Timman et al (1962). I f we restrain the cylinder from moving again, this be-comes
_ a4>;,ic
dn dn
Unlike the fluid domain our computational domain cannot be infinite, so we have to introduce artificial boundaries and proper boundary conditions. These conditions include the ra-diation condition mentioned earlier. A t the vertical bound-aries we choose this condition such that i t absorbs the out-going waves, but reflects the unwanted incoming waves. Because we have two outgoing waves on each side, i.e., we only consider T = {Uw/g) < 0.25, we get a second-order par-tial differenpar-tial equation
d d\ d a , - + C i — - - h C 2 — <|) = 0
dt dni \dt dn (10)
See also Romate (1989), who considers different angles of incidence i n 3D instead of different phase velocities. This condition can easily be adjusted for larger values of T . After substituting the phase velocities this equation reduces to
a% I — +U-±2
dt^ V T dndt dn" 0 (11)
On the right boundary the plus sign should be taken, on the left the minus sign. Here the steady potential (4) is used. We could also have used the double-body potential, but be-cause the artificial boundary is chosen relatively far away from the ship, these two are approximately equal.
At the artificial bottom of our domain we submit (j) to the condition
34)
dn = 0
This represents a real bottom, which would cause unwanted reflections. This is overcome by choosing the bottom rela-tively far away.
Movement of ship
To calculate the potential due to the movement of the ship, we have to find the forces acting on the ship. These first-order forces are given by
d<^
dt V(j) • V4) + V4)(X • VW^ nds
where ^ is the potential due to the incoming and the dif-fracted wave only, because the motion of the ship is yet un-known. The last term i n the integral is proportional to the movement of the ship and can be seen as a restoring coef-ficient i n the equations of motion. I f we do not consider in-coming waves but forced oscillation of the ship, we can cal-culate the added-mass and damping matrices, A and B , by fitting the first two terms of the force to the acceleration and the velocity: a^x F = - A — dt^ a x B — dt
The movement of the ship can then be calculated by solving the following differential equation:
a ' x a ' x a x
M— + A— + -B — +RX (12)
w i t h M the mass matrix and R the matrix of restoring coef-ficients. As was mentioned, the matrix R is partly built up by the third term of the first-order force, which is of second order i n the velocity. A t zero forward speed, only does not vanish.
These equations of motion can only be interpreted i n the frequency domain. The coefficients are frequency-dependent. Therefore, we restrict ourselves for the moment to harmonic functions.
Second-order forces
After we have calculated the movement of the ship, we can calculate the potential due to this motion. Then om- un-steady potential (6) is known, for both 4'i,ie sjid 4>d!/ v^ere al-ready known. Now we can calculate the averaged second-order forces, correct up to second-second-order i n the velocity, by the following formula:
(Fs) - p I ^ ( X - V ) ( ^ +v^-v<t>
dt
• (X • V)(V|(X • V)Vcf.) -h - Vt}) • V4) nds ) (13)
where 1,^ is the linearized relative wave height, which can be calculated by Bernoulli's equation and the movement of the hull.
Numerical algorithm
To solve the problem described above, we introduce the Green's function of the 2D Laplacian
G{x, ê ) = — I n \x-k\ (14)
By using Green's second theorem we have for the potential in the fluid domain
8(|)(x, t) with 3S (^(C, t) — ( X , I) - G(x, ^ 1 d r (15) dn dn. 8 = 1 x^S ~ xèdS 2 0 elsewhere
and dS the boundary of the domain.
For all the relevant quantities that can be obtained by the values of the potential on the boundary, only values of x on the boundary will be considered.
We discretize the integral equation (15) by dividing the boundary into intervals and assuming the quantities to be constant on such an interval. The collocation points are the midpoints of these intervals; see also Yeung (1973). This gives a set of linear equations which can be rewritten as:
D(|)(rt = E(|)„(Ö (16)
with D and E matrices built up by the Green's function and its derivative, and <j>n a vector containing the normal deriv-atives of (|). This set of equations could be solved i n time by iterating i n the following way:
• prescribe (j) or (j)„ on the boundaries at a time level t • solve (16) at this time level
• integrate the free surface and the artificial boundaries i n time to get new values for cf) or on the boundaries at the time level t + At
• go back to the second step.
Unfortunately this time integi'ation is unstable for eveiy size ofthe time step, although implicit schemes were used for the time integi-ation of the boundary conditions. This is due to the fact that after discretizing the space dependency, the free-surface condition contains several eigenvalues with real part gi-eater than zero for all non-zero velocities. For zero forward speed these eigenvalues become purely imaginary, so the above scheme can be used only i n that special case.
To achieve stable time integration, we developed a new method. I t is based on substituting the boundary conditions themselves i n (16) and integrating the resulting matrix equation i n time. Here we give the derivation i n the case that we make use of (8). I n the Appendix the formulas can be found for the general case when we take (7) instead of (8). First we make d^/dn explicit:
dn dxdt + u
,d%
dx' (17)
with subscripts denoting the time level. For convenience we split up all variables into a free-surface part and a remain-ing part:
.}>^=(4>j|ci>r)
By doing the same for the matrices and substituting the value of d^f/dn we get
The second-order time derivative can be discretized by a first-order difference
3 ^ dt'
1
m 2 (<i>i+i - 24,; + <t>f^i) + om
and the order derivative by its usual backwards first-order difference scheme. The spatial derivatives can all be approximated by second-order differences, while the values of the potential on the edges are extrapolated from the in-tersecting boundary. Applying this to the above equation and sorting the potentials according to their time-levels, we get
1
m
' E 2U
Af 0
4>,-+ E,4), The matrices Ci and Cg result from the discretization of the spatial derivatives. So now the unknown d4>f/dn has been eliminated. I n the derivation we used (4) while (5) is used in the Appendix, but the first one was used here i n order to keep the formulas short and understandable.
The same procedure can be repeated for the artificial boundaries. The second derivative i n the normal direction can be replaced by the second derivative i n the tangential direction by means of the Laplace equation. This leads to the following overall matrix equation:
Di4)i+i = Da.}», + Bs<^i-i + f;+i (18) with f a time-dependent vector resulting from the body
boundaiy condition and parts of the spatial derivatives. With the condition that the fluid is initially undisturbed, i.e., cj) = 0 and d^/dt = 0, this system can be solved.
Results
The numerical algorithm presented i n the previous sec-tion has been used to calculate the second-order forces and all other relevant quantities, such as added mass and damp-ing coefficients and movements of the cylinder. We studied the case that D = R, i.e., the draft is half the cylinder di-ameter, i n order to be able to compare w i t h literature. We will compare om- results with the results found by Zhao (1988). He was interested only i n small values of U and used a free-surface condition which is linear i n U.
We used a computational domain which had a length of three wavelengths on both sides of the cylinder and a depth of one wavelength. On the free surface an equidistantial gi'id was used. On the upstream side, 60 panels were used and on the downstream side more than 60 panels, because ofthe longer waves at this side combined w i t h the equidistantial gi-id. A t the bottom only 10 panels were used and on the artificial boundaries the interval size was twice the in1;erval size of the free surface. Larger elements on the artificial boundaries would give rise to diverging terms, because of the extrapolations needed at the corners of the domain. On the cylinder 40 panels were used to get enough accuracy for the pressure integration and the derivatives ofthe potential.
The derivatives of the potential on the cylinder were ob-tained by careful discretization. Most of the discretization was straightforward, but at the intersection of the cylinder with the free surface, the existence of stagnation points has to be taken into account. We used the value of the potential at this intersection to discretize the derivative. This value had to be extrapolated from the free surface. This is possible only i f the potential on the free surface is known accurately.
Fig. 3 Added mass in heave for six velocities Fig. 4 Added damping in heave for six velocities
Therefore the number of points on the free surface can only be reduced using a non-equidistantial grid.
The time integration was carried out over a time interval of 8 periods according to the frequency of encounter. On this interval 400 time steps were taken. The equations of motion were integrated using the implicit method of Crank-Nichol-son.
The accuracy of the results was influenced mostly by the number of panels on the free surface, especially for low fre-quencies. This can be explained by the fact that the com-putational domain becomes very large for these frequencies and the steady potential is nearly zero at almost every col-location point. By taking more panels, the region where this potential is of influence can be represented more accurately. The relative error can be estimated roughly at a few percent.
As the incoming potential we used cos(a)f — ^x)e* w i t h and w = W o + kU k = (Oo g
The following results were all calculated using the double-body potential (5), for the results with the steady potential (4) were unsatisfactory. This can be explained by the fact that the normal velocity on the ship hull should be zero for the steady potential, which is fulfilled by (5), but not by (4). We calculated the results for six different Froude num-bers: Fn = -0.035, 0, 0.035, 0.07, 0.105, 0.14, with the Froude number given by U/'VgR. The added-mass and damping coefficients agi-ee well with the results found by Zhao (1988), but for higher velocities the coefficients get influenced by this velocity (see Figs. 3-6). This is due mainly to the sec-ond-order free-surface condition. We see that the added mass i n heave is a symmetric function of U in Fn = 0, for the curves for Fn = ±0.035 coincide. The other coefficients seem to be asymmetric. Furthermore, we notice that the maxima i n the sway coefficients become more extreme when the ve-locity increases, which can also be concluded for the mini-mum in the added mass i n heave.
The added-mass and damping coefficients i n the case of zero forward speed were also checked w i t h Vugts (1968). The heave coefficients agree well, but the sway coefficients agi-ee only after we increase the number of panels on the free sur-face. This is because the large gradient occurring i n sway motion near the intersection of the hull and the free surface
can be represented accurately only by using very small panels. The problems could have been overcome by using a multi-sized grid, but this was not necessary to achieve our goal i n this investigation.
The coupling coefficients could not be checked, but they look rather convincing; see Figs. 7-10. Only the coupled damping coefficients remain of one sign; for the mass coef-ficients there seems to be a common zero for all velocities. Note also that the coupled coefficients become rather large for increasing Fn, so they cannot be neglected i n (12).
The results for the horizontal drift force are given i n Fig. 11. For low velocities they agi-ee well with the results Zhao (1988) found. For higher frequencies we see clearly that we included nonlinear effects i n U, for the increase i n the d r i f t forces lessens due to an increase i n the velocity. The extreme value of the d r i f t force is shifted to the left, for only the
fre-0.6 0.45 0.3 0.15 pd>-Fn = -.035 Fn = 0 Fn = .035 Fl = .070 Fn = .105 Fl = .14 Wo 0.375 0.75 1.125 1.5 V 3 Fig. 5 Added mass in sway for six velocities
0.375 0.75 1.125 Fig. 6 Added damping in sway for six velocities
-0.2 I
Fig. 7 Coupled added mass in heave for six velocities
-0.3
Fig. 10 Coupled added damping in sway for six velocities
quency of encounter matters for this value. To obtain these drift forces i t was very important to take second-order terms i n U into account, especially for larger U. I f one would not consider the large influence the velocity has on the added-mass and damping coefficients, the drift forces would be-come very extreme and unrealistic.
Figure 12 shows the results for the vertical second-order force. Here we see that the influence of the forward speed merely results i n a shift of the drift force, due to the fre-quency of encounter.
Unfortunately no measurements are available yet to ver-i f y our calculatver-ions.
the equations, which are correct up to second-order i n the current velocity U. This algorithm is stable for every veloc-i1:y and for every time step we tried. Therefore i t is a very powerful algorithm which gives hope for the existence of a similar algorithm i n 3D.
To test out the method we used harmonic functions. I f we consider a general time signal as an incoming wave, the equations of motion (12) have to be replaced by the proper ones as described by Ogilvie (1964). The effect of these equa-tions on our method will be the subject of future investiga-tions.
For low velocities the calculations agree well with the re-Conclusions
In this paper we tackled the time-dependent ship-wave problem. A numerical algorithm has been developed to solve
0.3 0.2 0.1 -0.1 \ \ Fn = -.03.') FJI = 0 Fn = .03.5 Fn = .070 Fn= .105 Fi= .14 Fn = -.03.') FJI = 0 Fn = .03.5 Fn = .070 Fn= .105 Fi= .14 Fn = -.03.') FJI = 0 Fn = .03.5 Fn = .070 Fn= .105 Fi= .14 Fn = -.03.') FJI = 0 Fn = .03.5 Fn = .070 Fn= .105 Fi= .14 0.375 •r.'l25 0.8 0.6 0.4 0.2 Fn = -.035 Fn = 0 Fl = .035 Fn = .070 Fn = .105 Fl = .14 0,375
Fig. 11 Drift force for six velocities
0.)
-0.1
-0.2
Fig. 8 Coupled added damping in heave for six velocities
1.125 Fl = -.035 Fn = 0 Fl = .035 Fn = .070 Fn= .105 Fl - .14 Fl = -.035 Fn = 0 Fl = .035 Fn = .070 Fn= .105 Fl - .14 Fl = -.035 Fn = 0 Fl = .035 Fn = .070 Fn= .105 Fl - .14 0.3 -0.3 -0.6 pgQ Fi = -.035 Fn = 0 Fn = .035 Fn = .070 Fn = .105 Fl = .14
suits Zhao (1988) found using the frequency domain and a linear free-surface condition.
References
M O R E T T I , G . 1964 Functions of a complex variable. PrenticeHall, E n -glewood Cliffs, N.J.
O G I L V I E , T . F . 1964 Recent progress towards the understanding and prediction of ship motions. Proceedings, Fifth Symposium on Naval Hydrodynamics, Bergen, Norway.
R A V E N , H . C . 1992 A practical nonlinear method for calculating ship wavemaking and wave resistance. Proceedings, 19th Symposium on Naval Hydrodynamics, Seoul, South Korea.
R O M A T E , J . E . 1989 The numerical simulation of nonlinear gravity waves in three dimensions using a higher order panel method. Ph.D. thesis. University of Twente, The Netherlands.
T I M M A N , R . A N D N E W M A N , J . N. 1962 The coupled damping
coeffi-cients of a symmetric ship. J O U R N A L O F S H I P R E S E A R C H , 5, 1, March. V u G T S , J . H . 1968 The hydrodynamic coefficients for swaying,
heav-ing and rollheav-ing cylinders in a free surface. Report No. 194, Delft Uni-versity of Technology, Delft, The Netherlands.
Y E U N G , R . W - C . 1973 A singularity-distribution method for free-sur-face flow problems with an oscillating body. Report No. NA 73-6, Col-lege of Engineering, University of California, Berkeley.
Z H A O , R . A N D F A L T I N S E N , 0. M . 1988 Interaction between waves and current on a two-dimensional body in the free surface. Applied Ocean
Research, 10, 87.
I I X
In this Appendix we will derive the matrix formulation of our time iteration as shown in the section "Numerical algorithm," but now for the free-surface condition (7). This condition was used to calculate the results shown before. The free-surface condition is
dxl dx^ ax' dt
fl> 1 B%
dxj J \dx^ gSzSi
By making dif/dn or dif/dz explicit we get
a'^\\ i [ i 2 — —- 1 3 — - H dx dxdt dx a . r dx
+
a4) = 0 atz = 0 dz dn 2g' \ \dxj I dndf} dx dxdt dx dx^ dx _^ ld^Vd% _^ a<|) 1 dx) dx' dx' dt 2 d% If we define lix) = -\U' dxand discretize the time derivative of d^/dn we get I \a<)> (gMY) dn I dn dn 1 (d'^ at}) d'<^ d^ a^4> 5<j) g 1 dt' dx dxdt dx dx' Bx ^ / a 4 ) y a ^ ( t > ^ a ^ f l < ) ) _ ^ a ^ \dx) dx' a.x-^ dt dx"
Again, we split up all variables into a free-surface part and a remaining part:
By doing the same for the matrices and substituting the value of a<j>//
dn we get a<j>l D^(t>f,,-+i + Ef I P(g^ïf dn 34) \ 1 + — + -i dn i-J pg L d'^ acj) a > dx dxdt dt a<j) a^4) a<j) /a<t,\ a"<t> d'i/d^ d'^ dx dx' dx \dx/ dx' dx' dt dx'
with
p = pix) = 1 -1
lix)
Now the terms between the square brackets can be discretized in the same way as was done in the example. So we get
Dr+-E, g [iMf m -1 0 C|>,-^i + Er^r.nMl + feAO-2 C C l ] {Mf (A«)J ErC(2iff,„,i - <|),-,„,i-l)
The matrix C arises from discretizing l/p(x) and is diagonal. The matrix C l is built up by the discretization of the terms with both a spatial and a time derivative. All the other derivatives are combined in the matrix Ca.
