Vorticity Fields due to Rolling Bodies in a Free Surface—Experiment and Theory
R.Yeung, C.Cermelli, S.W.Liao (University of California at Berkeley, USA)
ABSTRACT
Viscous effects are known to have significant influence on hydrodynamic forces on bluffshape bodies. Ocean structures in long waves and rolldamping arising from bilges of a ship hull are important examples. Some recent efforts to understand and model such realfluid effects in the presence of surface waves are described. As a canonical problem for studying ship roll, a surfacepiercing plate is oscillated sinusoidally in the free surface. A Digital Particle Image Velocimetry system (DPIV) was developed at the UC Berkeley ShipModel Facility to map out the flow details. In one regime, the flow exhibits an asymmetrical vortex shedding pattern. This steadystate pattern is observed to be systematically dependent on the initial conditions. DPIV and other experimental results were used to validate the predictive capabilities of the FreeSurface RandomVortex Method (Yeung & Vaidhyanathan, 1994). This gridfree Lagrangian method is briefly described. With some calibrations, the method reproduces very well the observed vortical structures, the measured forces and moments, and the wavegeneration characteristics. Application of the methodology to a rolling rectangular section with sharp bilges is also considered. The limitations and promise of the method are discussed.
1
INTRODUCTION
To be successful in predicting the motion of marine vehicles, particularly motion in the horizontal plane, sway, roll and yaw, one needs to be able to model viscous effects in a reasonably accurate way. In the last decade, much strides have taken place in the development of solution based on potentialflow theory that includes the effects of free surface (Yeung, 1982, Korsmeyer et al., 1988, Newman, 1992). While heave and pitch motions are generally considered to be well predicted, successful predictions of roll motion still rely on empirical estimation of roll damping, often on a trialanderror basis (e.g. Schmitke, 1978, Brown et al., 1983).
Progress in this area has been slow since the proper modeling of viscous effects would involve the solution of the NaviesStokes equations with a capability of modeling massive separation arising from the bilge corners or bilge keels. Further, one may also need to account for coupling coming from the freesurface.
A closely related problem involoving flowseparation effects occurs in offshore engineering where the bluff shapes of the structures lead to substantial “drag” and “inertia” forces in waves. Here, separation effects are commonly characterized by the “KeuleganCarpenter Number” KC= U_{m}T/D, where U_{m} is a characteristic velocity, T the time period of the flow, and D the body dimension. When KC is of O(10) or higher, realfluid effects are known to be significant, (Susbielles & Bratu, 1981, Sarpkaya & Isaacson, 1981). Typical applications have KC values ranging from O(10) to O(10^{2}).
The objective of the present paper is to report our recent efforts in developing a model for such separatedflow phenomena, particularly in the presence of surface waves. When a marine vehicle undergoes extreme roll motion, separation from the bilge corners, or from bilge keels are the dominant contribution to damping, whose magnitude critically determines the severity of motion, especially in conditions close to resonance. Thus any theoretical or computational model should be properly validated. This requires careful complementary experimental work which is also described in this paper.
It is worthy to mention that there are existing semiempirical methods for rolldamping estimation due to the hull surface and bilge keels (Himeno, 1981, Ikeda, 1977, Tanaka, 1960). Of course, one component of damping is surfacewave generation, which is normally well predicted by potentialflow theory. There have been efforts also to predict viscous damping without relying on empiricism, (Fink and Soh 1974, Brown and Patel 1985, Braathen and Faltinsen 1988, Cozens 1987, Downie et al. 1990). None has yet been able to model the interaction of hull geometry, vorticity generation and freesurface simultaneously.
Recently, Yeung & Vaidhyanathan (1994) have developed a FreeSurface RandomVortex Method (FSRVM) to model nearsurface flows for body geometry that may or may not have sharp edges. Their earlier works involved several important improvements (Yeung et al., 1993) over the Random Vortex Method (RVM) proposed originally by Chorin (1973). An efficient numerical scheme was devised to compute the mutual induction of discrete vortex elements. A boundaryintegral formulation was used to handle arbitrary body geometry. These new developments were extended to take into account of freesurface effects. The essence of this method is briefly described in Section 4. FSRVM does not require any prior knowledge of the location of the separation point. Since the method is gridfree, it can represent vortex structures of a rather wide range of scales. Accurate predictions at high Reynolds number is also possible, as illustrated in Section 5 by some direct comparisons of the numerical predictions with experimental results.
Although forces and moments are usually the primary quantities of interest, it is not always appropriate to establish the validity of a mathematical model on forces or moments alone. A more rigorous validation should involve comparing flow patterns also. At the University of California ShipModel Testing Facility, a Digital Image Particle Velocimetry (DPIV) was developed and implemented to provide quantitative results on the velocity and vorticity fields generated by a body moving in a free surface (Cermelli, 1995). Some of the major findings of such an investigation are reported in Section 3.
Because of the complexity of vortex patterns generated by a realistic hull model, a simplified “canonical” problem was introduced to study the roll motion of a ship section. The canonical problem chosen is that of a surfacepiercing plate rolling periodically about the free surface. The flow is primarily two dimensional so that comparison with applicable theoretical results is immediately possible. For the typical frequency range of interest, this rolling plate is designed so that surfacewave generation and and vorticity generation are important. Naturally, besides providing a validation on our theoretical model, this canonical experiment is closely relevant to the understanding of the performance of bilge keels, the prevalent passive device for reducing roll motion. In the last section of the paper, we present predictions of added inertia and “equivalent” linear damping for a rectangular hull section and compare them with some existing data.
2
THE ROLLINGPLATE EXPERIMENT
2.1
Design of experiment
A special model was designed and fabricated to study the vortical structures generated by a plate undergoing rolling motion in a free surface. This experiment was carried out at the ShipModel Testing Facility of the University of California at Berkeley. The towing tank is 61m(L)×2.44m(W) ×1.52m(D). The plate is oscillated by a hydraulic piston. Flow visualization can be conducted through four observation windows along tank walls.
The plate was built out of oneinch thick acrylic and stiffened longitudinally. It has a natural frequency of 4 Hz, while the range of forcedmotion frequency being investigated is less than 0.6 Hz. The clearance between the plate vertical edges and the tank walls was kept to a minimum (0.64cm) to reduce end effects. The plate was hinged at the free surface. Its draft was primarily set to 30.48cm (12 in). Flow visualizations were conducted also with a 15.24cm (6 in) draft.
A plate with a tip that was rectangular in section was found to yield force and moment measurements that were not repeatable, primarily as a result of chaotic vortex shedding from the bottom edges of the plate. The tip of the plate was then modified by the addition of a half circle rod at the bottom of the plate (see Fig.1). This greatly improved the consistency and repeatability of force records.
Plate rotation was achieved by applying horizontal motion to a rod which was hinged to the plate on one side and to (a random motion) hydraulic piston on the other side (Fig.1). A software controller was developed to ensure speedy
response of the piston to any prescribed function of time. The software controls the plate angle α(t) according to the following prescribed time function:
(1)
where ω=2π/T is the angular frequency of the prescribed oscillation and α_{0} the roll amplitude. α_{0} can be either positive or negative, depending on whether the plate starts with a right or left swing. The “smoothness parameter”, s, determines the number of periods needed to reach a sinusoidal steady state of α(t). The plate motion always starts with the vertical position, α=0.
For later purposes, it is helpful to introduce a phase angle Φ to define the plate position. Φ equals to 0° is defined to be when the plate is vertical (α=0) and the plate tip is moving toward the right (i.e., increasing α). Φ=90° when the plate reaches its maximum angle of swing (α=α_{0}). Φ=180° represents the instant of α=0 when the plate tip is swinging toward the left (decreasing α). Thus, Φ=270° corresponds to the plate at its minimum angle position (α=–α_{0}).
Force transducers were installed at the two side hinges holding the plate and at the hinge connected to the driving rod to measure the moment and the horizontal and vertical forces. The equations of motion of the plate were developed. The desired hydrodynamic forces and moment can be obtained by deducting contributions of the hydrostatic component (computed from the plate roll angle) and of the inertia components from the measured values. Output from the force transducers and from a wave gage were digitized at 50 Hz and recorded on a 486PC. In addition,
the pistonlocation feedback signal and the wave elevations at two points were recorded.
2.2
FlowImage Capturing and Analysis
The Digital Particle Image Velocimetry (DPIV) System set up for this experiment is shown in Fig.2. The light source consists of a 10W argon laser generates a 570 nm wavelength (green) beam of approximately 2 mm in diameter. To avoid light refraction at the free surface, fiberoptics are used to direct the laser light into the tank. The beam emerges into a cylindrical lens to form a divergent light sheet shining on the rolling plate from the underside.
Two video systems are used to record successive images of the flow. One is composed of a blackandwhite CCD camera with RS170 output, and a Sony laser disk video recorder. The second system is a Hi8 Canon Camcorder that records images on a video tape with high fidelity. A speciallydesigned electromechanical shutter is installed between the laser head and the fiberoptic coupler to control the exposure time of the cameras.
For flowfield seeding, lowcost Pliolite particles of the Goodyear Company are used. It has a specific gravity of 1.05. Pliolite granules were grind down to particles of size about 300 µm. With some presoaking, they are less buoyant, and their reflective properties make them quite
suitable for DPIV.
A threestep procedure was developed to compute the velocity and vorticity fields from pairs of successive digitized images at 1/30sec apart. In the first step, a wellknown crosscorrelation algorithm of Willert and Gharib (1991) is implemented to yield a good estimate of the flow field. In the second step, spurious vectors are discarded or replaced by reinterpolated neighboring values. In a third step, suggested by Fabris (1994), a modified version of the crosscorrelation algorithm is performed to obtain a more accurate representation of the flow field. In theory, these last two steps can be iterated until a converged velocity field is obtained. In practice, results were found to be satisfactory after only one iteration.
The algorithm implemented in the first step consists of the following:

Select an interrogation window (F_{1}(m,n) in digital form) on the first image where m and n are pixel indices in the horizontal and vertical directions and F_{1}(m,n) is the grey scale value at (m,n).

Select an interrogation window (F_{2}(m,n) in digital form) on the second image taken a time Δt later. The two interrogation windows (F_{1}(m,n) and F_{2}(m,n)) correspond to the same location on each image.

Compute and Fourier transforms of F_{1}(m,n) and F_{2}(m,n) respectively, using FFT.

Multiply by (complex conjugate of to obtain the function Φ(u,v) and compute its inverse Fourier transform. Illustrative representations of the pixel functions F_{1}(m,n), F_{2}(m,n) and (m,n) are shown in Fig.3.

The exact location of the highest peak of ((m,n)) is obtained with subpixel accuracy by fitting a parabola around the maximum of ((m,n)). The shift in horizontal and vertical direction (d_{x},d_{y}) is the distance from the peak to the origin. Horizontal and vertical velocity are obtained by dividing (d_{x} and d_{y}) by the time step between the two frames.
The process above is repeated after shifting the location of the interrogation window horizontally or vertically by a small amount (16 pixels). One vector corresponding to the velocity at the center of each interrogation window is obtained every 16 pixels. Because of the large size of the interrogation window (64×64 pixels), the velocity is usually strongly smoothed out. The advantage of using a large interrogation window is that higher velocities can be resolved better. In high velocity gradient regions, spurious vectors are sometimes obtained; they are removed in the second step described next.
In high velocity regions (near the vortex cores), spurious data is removed manually and velocity is estimated based on the streak length (usually 5 to 10 pixels). Accuracy of the results is not critical at this stage. Only a “reasonable” estimate of the velocity field is required since “fine tuning” is carried out in the third step. A graphic interface showing simultaneously the digitized images and corresponding velocity field was helpful in speeding up this “manual” process. Different algorithms based on statistical removal and reinterpolation of spurious data were not satisfac
tory due to the high density of erroneous data in some regions.
In the third step, smaller interrogation windows (16×16 pixels) are used and the interrogation window on the second image is shifted by an amount corresponding to the shift (d_{x},d_{y}) determined by the previous steps. Crosscorrelation of the two interrogation windows is then carried out using a FFT based algorithm similar to that followed in step 1. This technique allows high velocities be resolved without sacrificing spatial resolution.
A thorough error analysis (Cermelli, 1995) has shown that peak vorticity may be in error up to 20 % near the vortex core and less than 10 % in other flow regions. Figure 4 shows a digitized image of a typical data frame.
3
VORTICITY STRUCTURES— EXPERIMENTAL RESULTS
For the specific setup described in Section 2, the “Rolling Plate Experiment” can be characterized by the three dimensionless parameters. Let L be the draft of the plate. An appropriate Reynolds number is based on the maximum tip velocity of the plate and the distance traveled by the plate tip:
(2)
where v is kinematic viscosity.
The KeuleganCarpenter number (KC), a traditional parameter in the study can be defined based on the plate thickness t_{P}: :
(3)
The third parameter is the nondimensional frequency, which is defined by:
(4)
Limited by the capacity of the force transducers, the range of is confined to [.28, .63] for the experiment. Even so, this is a range that is quite close to typical rollresonance of ships.
For most runs, L was set at 31.75cm (12.5 in) and the following parameters were varied: the roll amplitude, α_{0}, (from 5° to 15°), the roll period, T, (from 1.8 to 4 seconds), the number of periods for the plate angle to reach full stroke s and the direction of the first swing of the plate, sign of α_{0}. Additional runs were made with L set at 16.51cm (6.5 in) and with the same motion parameters. Flow visualizations were carried out also at different stations across the tank to study the wall influence. The flow is partially turbulent in all runs, and many small threedimensional vortical structures interact in an intricate manner. However, after several periods of oscillation, a harmonic steady state can be observed for the main vortical structures.
3.1
Asymmetrical flow regime
A significant flowfeature that was uncovered during the experiment is somewhat counterintuitive and will be referred to henceforth as the “asymmetrical flow regime”. Figure 5 shows results that correspond to the case of KC=20.6 (L= 31.75cm), with a starting swing to the left. In this and other subsequent plots, the following conventions are adopted. The x and y scales represent linear dimensions in mm. A vector of 20 mm length in this scale represents a velocity magnitude of 170 mm/sec. Vorticity contours are shown at increment of 3 vorticity units (3 sec^{–1}), with solid lines taken as counterclockwise and dashed lines as clockwise.
Referring to Fig.5, we note that during the first half of a typical period when the plate tip moves toward the left, with a decreasing, a counterclockwise vortex rolls up in the wake of the plate (Fig.5a, Φ=158°), which is to be expected. When the plate stops at the minimum roll angle (–α_{0}), a small jet is created between the backside of the plate and the existing vortex, generating as a consequence a clockwise vortex (Fig.5b, Φ=280°). When the back swing takes place (α increasing), the strength of the second vortex increases quickly. When the two vortices
are of similar strengths, they move away from the plate as a vortex pair at about 45° angle with respect to the vertical (Fig.5c, Φ=360°). As the vortex pair moves away, it breaks down into turbulent motion and the vorticity intensity in the light sheet plane decreases rapidly. The vortex movement of the return half cycle is thus distinctively different from that associated with the forward half cycle. Because of this asymmetry, the magnitude of positive and negative peaks of the moment and horizontal force are noticeably different, as illustrated by Fig.6. In this figure, the sway force time history F_{x} is shown together with α(t). The periodic asymmetric behavior is evident. This asymmetrical flow regime was observed for KC in the range [13.7, 20.6] and R_{e} in the range [9,600, 28,900]. They correspond to all runs with L=31.75cm and α_{0} equals to 10° or 15° and to values of between .28 and .63.
For flows without a free surface, a somewhat similar phenomenon was reported by Singh (1979) and more recently by Sarpkaya and O'Keefe (1995). In the present work, it was possible to identify the factor that determines the direction toward which the vortex pairs are shed. It was found that it depends almost entirely on the direction of the first swing. Out of 32 runs resulting in vortex pairs moving toward the left, 31 started with a left stroke. Out of 14 runs resulting in vortex pairs moving toward the right, 13 started with a right stroke. This excellent correlation between direction of shedding and direction of the first stroke is related to the use of a rounded bottom edge (with a rectangular edge, the direction of shedding was random). One may surmise that a plate edge with a triangular tip would behave similar to the roundedbottom case.
Vorticitycontour plots of the transient state have shown that a vorticity imbalance originating during the first few strokes of the plate is responsible for determining the direction toward which the vortex pairs move after establishment of a steady state. In some cases, this can take as many as ten periods of motion. We have observed a very nonlinear phenomenon, in which the longtime behavior of the system depends on very small perturbations generated during the first few periods of motion—a characteristics of “bifurcation”.
In most runs, the light sheet was located at the tank centerplane (half way between the walls). Several experimental runs were made with the light sheet located at different transverse stations across the tank. Figure 7 shows the flow two feet from the tank wall for a run with T=2 seconds and α_{o}=10°. The flow observed is similar to the flow at the tank centerplane. A vortex pair is shed toward one side of the plate. The direction of shedding also corresponds to the direction in which the plate tip moves at the first stroke. Flow visualization taken at a distance of six inches from the tank wall do show strong “endwall effects”. The wall inhibits the vortexpair from moving away. Instead, they hover around the tip of the plate.
3.2
Symmetrical flow regime
For smaller roll amplitude (α_{o}=5 deg) and for all the runs with a 16.51cm draft, a different flow regime is obtained. It is more consistent with one's intuitive expectation (see Yeung & Cermelli, 1996).
In a typical scenario, a vortex is generated on the plate back side (say the right side as the plate tip moves toward the left). This vortex remains attached to the plate and as the plate reverses direction, a smaller vortex of opposite sign is generated on the left side of the plate. The resulting vortex pair is not strong enough to move away from the plate, and the growing vortex eventually absorbs the earlier one as the plate's backswing continues.
In all, the flow patterns corresponding to motion of the plate in one direction are mirror images (about the vertical axis) of those corresponding to motion of the plate in the other direction. In terms of the nondimensional parameters, this “symmetrical flow regime” corresponds to KC between 6.9 and 10.0 and R_{e} between 3,200 and 9,600, with varying from .28 to .58.
4
THEORETICAL MODEL
A very powerful method for modeling viscous flow in the presence of surface waves is the FreeSurface RandomVortex Method (FSRVM), developed recently by Yeung and Vaidhyanathan (1994). The RVM is a gridfree formulation which can capture vortex scales of a wide range. Validation on this model has been carried out for a number of unbounded fluid problems (Yeung et al., 1993). The experimental results of the last section provide an excellent opportunity to validate and calibrate the theoretical model. Towards this goal, the FSRVM code was improved and adapted to solve the rolling plate problem. The original formulation of this method can treat bodies of arbitrary shapes and three degrees of freedom of motions of an arbitrary type. Hence, the theoretical development outlined below is applicable to problems of bodyfreesurface interaction of a general nature.
4.1
Viscous flows with a free surface
The fluid domain in question is defined in Fig. 1, with the boundary consisting of the free surface the body surface and the open boundary the latter being an artificial boundary introduced to truncate the computational domain. The Random Vortex method solves the NavierStokes equations using a vorticitystream function formulation. Thus for an incompressible Newtonian fluid, the governing equation for the vorticity vector ξk=∇×u is the vorticity transport equation:
(5)
where D/Dt is the material derivative. The stream function satisfy the Poisson equation:
∇^{2}ψ=–ξ. (6)
For material points on the free surface the kinematic freesurface boundary condition can be written in the following form:
(7)
where x is the location of a material point and u is the velocity at that point. The dynamic boundary condition states that the stress should vanish at the interface (see Wehausen & Laitone, 1960). However our previous experiences with full NavierStokes solution in the presence a free surface (Yeung & Ananthakrishnan, 1992a,b) indicate that the shear layer due to the fullstress condition on the free surface is normally of secondary importance (compared with bodygenerated vorticity). Thus an effective dynamic boundary condition is simply an inviscid one:
p=0, (8)
which states that the gauge pressure vanishes.
On the body the no leakage and no slip boundary conditions apply:
u=v_{b}, (9)
where v_{b} is the velocity at a boundary point x on At any parametric point s along the body contour, v_{b} is given by:
v_{b}(s)=U_{b}+Ωk×(x(s)–x_{0}(s)), (10)
where U_{b} represents translation of a reference point (x_{0}(t),y_{0}(t)) on the body and Ω(t) is the angular velocity about x_{0}.
At the open boundary at a sufficiently large distance from the body, the fluid motion is assumed to be minimal. Hence the pressure can be taken as hydrostatic,
p=–ρgy. (11)
4.2
Vorticity representation by blobs
The twodimensional vorticity field of Section 4.1 is represented by a collection of regions of concentrated vorticity called “blobs”, i.e.
(12)
(13)
Here x is an Eulerian field point in the fluid domain or on the boundary at time t. Γ_{i} is the circulation associated with each vortex blob (of which there are N) centered at x_{i}. The function G describes the distribution of vorticity within each blob whereas F is related to G by integration. The ‘Chorin blob', which consists of a step distribution of vorticity in a finite core radius, was successfully employed in our earlier works (Yeung et al., 1993). Outside of the blob core, the flow is irrotational.
To obtain the complete velocity field, we let Ψ in be decomposed into two components:
Ψ=Ψ_{v}+Ψ_{h}, (14)
It follows from the Poisson equation (6) that
∇^{2}ψ_{v}=–ξ, (15)
∇^{2}ψ_{h}=0. (16)
ψ_{v}, the rotational component, is given by (Eqn. 13), which is known everywhere once the blob strengths Γ_{i} are known. On the other hand, the irrotational component ψ_{h} has to be solved to satisfy Laplace's equation, subject to the various boundary conditions on This nontrivial hydrodynamic problem is discussed in Section 4.4.
4.3
Vorticity diffusion and convection
Per Chorin's (1973) original algorithm, the vorticity transport equation Eqn. (5) can be split into two equations, each more tractable by itself than the combined equation. This is a fractional step method. At the first and second half of each time step, the convection part and the diffusion part of the equation
(17)
(18)
are solved successively.
Eqn. (17) corresponds to the convection equation of an inviscid fluid, with the blobs taken as material (lagrangian) points, i.e.,
(19)
where the velocity field u has contributions coming from both the freesurface and the body motion. The free surface effects can be included in this convection step.
The diffusion equation (Eqn. 18) can be simulated by the random walk algorithm (Chorin, 1973 and Einstein, 1956). In this algorithm, each of the blobs is given an independent random walk in both the x and y directions with standard deviation being the time step over which Eqn. (18) is integrated. The numerical implementation of this ‘integration' is simply:
x_{i}(t+Δt)=x_{i}(t)+(η_{x},η_{y}), (20)
where η_{x}, η_{y} are chosen from a set of random numbers with Gaussian distribution whose mean is zero and standard deviation The errors introduced by this method are statistical in nature and decrease with increasing number of blobs. The effect of the free surface is neglected in this diffusion model.
In many such computations, longtime solutions are often necessary in order to capture the flow features of interest. Equation (19) represents the heaviest computational burden since by the BiotSavart law, each time step would involve N^{2} interactive computations among the blobs. As N increases, the solution process slows down very quickly. Fortunately, by using the fast multipole method of Carrier et al. (1988), we were able to implement an O(N) interactive algorithm which greatly enhanced computation efficiency (see Yeung & Vaidhyanathan, 1993).
4.4
Solution of the stream function ψ_{h}
ψ_{h} satisfies the Laplace equation, which is most conveniently solved by a boundaryintegral equation technique. For free surface flow problems, the following complexvariable formulation following the works of Vinje and Brevig (1980) and Grosenbaugh and Yeung (1989) is especially attractive. Introducing the conjugate function of ψ, i.e., a velocity potential , we can write the timedependent complex potential β as:
(21)
Here z_{i} and Γ_{i}, i=1,.., N, are the complex coordinates and circulations of the N vortex blobs present in the fluid domain. For any field point outside of each individual blob core, the complex velocity can be written similarly as:
(22)
At any given time t, the location of the body boundary is prescribed while the location of the free surface boundary is determined from its previous location by using the kinematic boundary condition (Eqn. 7). β_{h} is sought as an analytic function of z, subject to the inviscid freesurface condition Eqn.(8). In terms of , this condition can be written as:
(23)
The first term on the right hand side may be derived from Eqn. (21):
(24)
where, denotes the real part of the expression. In Eqn. (24), ż is the complex velocity at a Lagrangian point z located at the free surface and is the complex velocity of the ith vortex. Thus, Eqn. (24) can be used to compute the value of _{h} at the new location of the free surface, which is determined by a time integration of Eqn. (7). A Dirichlet condition for _{h} is therefore available at every discrete time instant on
On the open boundary, we approximate Eqn. (11) by
(25)
which can be expanded as:
(26)
Using fixed (Eulerian) node points on the open boundary, we can integrate Eqn. (26) to obtain a Dirichlet condition for _{h}.
Finally, the normalcomponent of the boundary condition on the body (Eqn. 9) can be implemented as a specification of ψ_{h} on . By (10), it is simple to show that:
(27)
where K(t) is an arbitrary timedependent constant and R_{0}=x–x_{o}. Note that the value of ψ_{v} on the body can be evaluated from Eqn. (13).
In summary, we have obtained a boundary value problem for β_{h}, where _{h} is specified on a portion of the boundary contour while ψ_{h} is specified on the remaining portion. This allows us to solve for either ψ_{h} or _{h} on all the boundary points using a boundary integralequation method based on Cauchy's integral formula:
(28)
Details on the treatment of such an integral equation can be found in (Grosenbaugh & Yeung, 1989). Once β_{h} is determined on the boundary, the velocity w_{h} at any interior point can be computed by using the following Cauchy 's integral:
(29)
So far, we have mentioned that the vorticity field is represented by blobs, without explaining how they are generated. In the RVM, the blobs are introduced to the flow domain through the noslip boundary condition. At the end of every convection step (i.e., the integration of Eqn. (5)), the Poisson equation (6) is solved, and then the fluid velocity is computed along the body surface . The tangential component of this velocity will not satisfy the noslip boundary condition (Eqn. 9). The boundary condition can be satisfied by generating vortex blobs of appropriate strength (circulation) and radius at the body panels (Chorin, 1973). It is then released into the flow during the solution of the diffusion equation (18).
4.5
Hydrodynamic pressure and forces
When the vorticity field is approximated by point vortices, inviscidflow results can be used to compute the forces and moment on the body (Yeung et al., 1993). In terms of the complex velocity potential, the pressure on the body is given by:
(30)
where ∂_{v}/∂t is given by the last equality of Eqn. (26). An integral equation similar to Eqn. (28) is used to solve for ∂_{h}/∂t (Vaidhyanathan, 1993). Forces and moments on the body can then be obtained by direct integrations of the pressure. This model does not compute the tangential stress on the body surface, which are typically insignificant if the flow has substantial separation.
5
THEORY VS. EXPERIMENT—COMPARISON OF RESULTS
The fully nonlinear FSRVM algorithms in the last section was implemented with two modifications. First, it was found that a linearized, instead of the fully nonlinear, freesurface condition was adequate for the range of parameters covered during the experiment. Specifically, Eqn. (7) and Eqn. (23) are now replaced by
(31)
For longtime simulations, this reduced computational time to about one half of that for the nonlinear case.
Second, a decaying model for the blob strengths was implemented so as to bound the number of blobs to O(10^{5}). The evolution in time of blob strength was chosen to be the following function:
Γ_{i}(τ)=Γ_{i}(0)[1–tanh(3τ–6)]/2, (32)
where τ=t/T is the time (in number of periods) after the introduction into the flow of blob i with intensity Γ_{i}(0). If the ratio Γ_{i}(τ)/Γ_{i}(0) is less than 10^{–3}, the blob is removed from the flow. With this decaying function, the strength of each blob remains almost unchanged for 1.5 period and then decays rapidly during the following period. Implementation of a decaying model results in substantial savings in computational efforts. Typically, 6 to 8 periods of motion are needed to obtain a periodic steady state. Further, this model can reflect in part the decay of vorticity due to turbulent breakup. The decaying function was calibrated by observing visually that vortices typically break down on the average of about 2 periods. Note that this model is not meant to be a turbulence model.
Several runs were made with all parameters set to those of the actual experiment. Flow structures obtained numerically turned out to behave somewhat differently from those observed experimentally. A downward adjustment of the Reynolds number greatly improved the match between numerical and experimental results.
Two hypothesis may explain this discrepancy. First, the number of blobs could be insufficient. This limitation arises from computer resources. A finer spatial and time discretization increases the number of blobs in the flow field. Flow configurations based on a larger number of blobs and the actual Reynolds number corresponding to the experiment showed much improved match. However because of staggering computational demand, simulations could not be run long enough to achieve a steady state.
The second hypothesis is that turbulence, which is not taken into account by the twodimensional FSRVM model, accounts for discrepancy between prediction and observed flow features. Eddy viscosity, in the actual experiment, may contribute to an increase in the “apparent ” viscosity of the fluid. A quantitative estimate of the relative influence of each of these factors is difficult to obtain. Some ongoing research aimed at implementing turbulence models on parallel computers may yield a better understanding.
A comparison of the experimental and numerical vorticity fields is shown in Fig.8 for the asymmetrical flow regime. This corresponds to case of R_{e}=19, 200, KC=13.7, α_{0}=10°, =0.565, with the numerical method run at R_{e}=6, 400. Evidently, very good agreement is observed, reproducing all interesting features of the flow discussed in Section 3.
Force and moment comparisons for this run are shown in Fig.9. The top figure of Fig.9 shows that the vertical (heave) force exhibits a double harmonic behavior. This peculiar behavior was first observed in some earlier computations by FSRVM and was later confirmed by experiment. The behavior is attributable to the nonlinear body boundary conditions, i.e., Eqn. (9) is satisfied at the instantaneous position. If the fluid were unbounded, it can be shown that the horizontal and vertical hydrodynamic force (F_{1}, F_{2}) and the hydrodynamic roll moment M_{6} are given “exactly” by
(33)
where the µ_{ij}'s are the hydrodynamic added masses of the plate. Since the coupling hydrodynamic coefficient µ_{16} does not vanish for a plate hinged at the end, the doublefrequency behavior of F_{2} is dominant. This behavior could not be predicted within the context of linearized motion theory which would include only linear terms in velocity and acceleration.
Figure 9 also shows the hydrodynamic roll moment for the same case. The agreement is very encouraging, representing an important improvement over potentialflow results. The wave elevation obtained at two locations, x=±355 mm, are compared with the predicted results in Fig.10. Excellent agreement is observed. Since the wave amplitude is directly related to wave damping, it provides a validation that the theoretical model can capture both viscous and wave effects well.
Figure 11 displays the details of the vortical structures generated under the free surface by the plate and the freesurface elevation for the case of α_{o}=–10°. Each dot represents a blob from the FSRVM solution. These six “snapshots” bracket about one cycle of oscillation.
A limitation of the FSRVM code is that the symmetric and attached vorticalflow features in the regime described by Section 3.2 could not be reproduced, even though the force and moment predictions were in good agreement with measured values.
6
HYDRODYNAMIC COEFFICIENTS OF ROLLING CYLINDERS
As an illustration of the application of the theoretical model to study ship sections of a general shape, we consider the problem of a rectangular section with relatively sharp bilges which undergoes prescribed rolling motion of the form:
α(t)=α_{o} sin(ωt) (34)
The specific geometry has a beamtodraft ratio B/D of 2.0. This shape was used in the experiments of Vugts (1968) which has a physical beam of 0.4m and a rounded bilge radius of 2.5mm. The frequency parameter is defined as before (Eqn. 4). The experimental Re number is 49,000. Guided by the rolling plate results, we adjust Re downwards to 25,000 in the numerical simulation.
The solution of this nonlinear problem at t/T= 5.75 for =0.5 is shown in Fig.12. The complexity of vortex array shed from the corners after only several periods is evident from this plot. The solution method is sufficiently robust that as many as eight periods of motion have been successfully run to allow one to examine the “steadystate” behavior. Figure 13 shows the time history of the heave force and roll moment due to this prescribed α(t). Also shown is the roll angle itself. Note the doubleharmonic behavior of the force as in the case of the rolling plate. There exists also a net “setdown” to pull the section downwards. The continuous shedding of vortices and interaction with existing vortices produce a roll moment signature that is quite complex. Even so, the overall pattern is still reasonably periodic.
If one adopts an equivalent linearized characterization of the roll moment as Vugts,
(35)
where µ_{66} is the added moment of inertia, λ_{66} is the total damping, a Fourier analysis of the moment signature will yield the µ_{66} and λ_{66}. The “moving window” analysis based on a window of one period yields rather steady values of added inertia and damping coefficients (see Fig.14). It is possible to separate the “equivalent” damping to linear and quadratic components (Yeung & Cermelli, 1996), but this would not be pursued here.
Simulations were carried out for this particular hull section so that the equivalent added inertia and damping could be extracted over a frequency range. The roll amplitude α_{o} is taken to be a
relatively small value of 0.05 radians (2.86 degrees). The range of considered was [0.6, 1.25]. The nondimensional added moment of inertia and damping are shown in Figs.15 and 16 respectively. In these figures, and are nondimensionalized according to and where ∇ denotes the sectional area of the hull.
A comparison of these results with existing data of Vugts shows that the inertia coefficients are overpredicted. However, the inertia coefficients are consistent with inviscid theory. In unbounded flow, our experience indicates that added masses are not strongly affected by flow separation. The experimental results of Vugts suggest otherwise. More extensive experimentations are needed to resolve this issue. The allimportant quantities of wave and viscous damping are well predicted. In particular, the deficiency of inviscid theory is evident. The theoretical model thus holds much promise.
7
CONCLUSIONS
The flow generated by a body undergoing periodic rolling motion was studied experimentally and theoretically. A “rolling plate in a free surface” is taken as the canonical problem. A Digital Particle Image Velocimetry (DPIV) system was developed and used to obtain accurate velocity and vorticity maps. For the range of parameters studied, when KC is larger than 13, it was found that an asymmetrical periodic steady state is reached. In this state, a vortex pair is shed toward a “preferred side” of the plate in one half
of each period of swing. This directional preference was identified to correlate with the initial direction of the swing and could be caused by the gradual increase of motion amplitude at the start of the swing. This graduality generates vortex pairs of unequal strengths. For smaller KC number, the flow patterns exhibit features that are symmetrical with respect to the vertical axis during successive half periods, a more intuitive result.
The theoretical model used is one based on solving the NavierStokes equations with an inviscid freesurface condition. The FreeSurface RandomVortex Method (FSRVM, Yeung & Vaidhyanathan, 1994) uses a gridfree EulerianLagrangian formulation. It is capable of capturing vortical structures of a wide range of scales. The method is not particularly effective in predicting features observed in the symmetrical flow regime. For the asymmetrical flow regime, i.e., large KCnumber flow, the FSRVM provides accurate predictions. Even so, in order to obtain excellent matching, the Reynolds number used in FSRVM simulations had to be adjusted downwards. It is surmised that this was due either to an insufficient number of blobs in the calculations or to a lack of turbulence modeling in the method. Neither represent obstacles that are insurmountable to overcome in the future. The FSRVM method also yields satisfactory force and moment predictions. In particular, a demonstration is made to compute the viscous damping of a rectangular hull section in roll. The predicted damping agree well with some existing laboratoryscale results, a considerably improvement over computations based on inviscidfluid theory.
It is also worthy to note that the heave restraining force generally appears as a double harmonic of the rollmotion frequency, a consequence of the nonlinear body boundary condition, not viscous effects. This characteristic is quite prominent in larger angles of roll, experiment and theory.
The extensive experimental and numerical works reported here have been complementary to each other. A continuing understanding of this subject will require such parallel efforts.
Acknowledgements
We are pleased to acknowledge the support of the Office of Naval Research under Grant N00014– 91J1614 & N00014–95–1–0980, with Jim Fein as the program officer. Partial support from grants of the University of California, the Lawrence Livermore Laboratory, and the Shell Foundation are also gratefully acknowledged. Thanks also go to Dr. M.Vaidhyanathan (3Es), Prof. P.Ananthakrishnan (FAU), Profs. O.Sãvas and D.Liepmann (UCB), who had assisted in many valuable ways during the course of this research.
References
[1] Brown, D.T. EatockTaylor, R. & Patel, M. H. ( 1983). “Barge motion in random seas—a comparison of theory and experiment”, J. Fluid Mech, 129, pp. 385–407.
[2] Brown, D.T. and Patel, M.H. ( 1985). “A theory for vortex shedding from the keels of marine vehicles”, J. Engrg. Math., 19, pp. 265–295.
[3] Braathen, A. and Faltinsen O.M. ( 1988). “Application of a vortex tracking method to roll damping”, Advances in Underwater Technology, Ocean Science and Offshore Engineering, 15, Society of Underwater Technology. pp. 177–193.
[4] Carrier, J.; Greengard, L. and Rokhlin, V. ( 1988). “A fast adaptive multipole algorithm for particle simulations”, SIAM J. Sci. Stat. Comput., 9, pp. 669–686.
[5] Cermelli, C.A.. ( 1995)., “Vortical flows generated by a plate rolling in a free surface”, Ph.D. dissertation, Dept. of Naval Arch. and Offshore Engrg., Univ. of Calif., Berkeley.
[6] Chorin, A.J. ( 1973). “Numerical study of slightly viscous flow”, J. Fluid Mech., 57, pp. 785–796.
[7] Cozens, P.D. ( 1987). “Numerical modelling of the roll damping of ships due to vortex shedding ”, Ph.D. dissertation, Dept. of Aeronautics, Univ. of London, UK.
[8] Downie, M.; Graham, J. and Zheng, X. ( 1990). “Effect of viscous damping on the response of floating bodies,” Proc. 18th Symp. on Naval Hydrodyn., Ann Arbor, Michigan, pp. 149–155.
[9] Einstein, A. ( 1956). Investigation on the theory of Brownian movement, Dover, New York.
[10] Fabris, D. ( 1994). Private communications.
[11] Fink, P.T. and Soh, W.K. ( 1974). “Calculation of vortex sheets in unsteady flow and applications in ship hydrodynamics,” Proc. 10th Symp. on Naval Hydrodyn., Cambridge, MA, pp. 463–491.
[12] Grosenbaugh, M.A. and Yeung, R.W. ( 1989). “Nonlinear freesurface flow at a twodimensional bow”, J. Fluid Mech., 209, pp. 57–75.
[13] Himeno, Y. ( 1981). “Prediction of ship roll damping—state of the art” Dept. Naval Arch. & Mar. Engrg., Univ. of Michigan, Rep. no 239.
[14] Ikeda, Y., Himeno, Y., and Tanaka, N. ( 1977). “On eddymaking component of roll damping force”, J. Soc. Naval Arch. Japan, 142, pp. 54–64 (in Japanese).
[15] Korsmeyer, F.T., Lee, C.H., Nemwan, J.N., & Sclavounos, P.D., ( 1988) “The analysis of wave effects on tensionleg platforms” Proc. 7th Int'l. Conf. of Offshore Mechanics & Arctic Engineering (OMAE), Houston, TX.
[16] Newman, J.N. ( 1992). “Panel methods in marine hydrodynamics”, Proc. 11th Aust. Fluid Mech. Conf., Tasmania, Australia.
[17] Sarpkaya, T. and Isaacson, M. ( 1981). Mechanics of wave forces on offshore structures, van Nostrand Reinhold Company.
[18] Sarpkaya, T. and O'Keefe J.L. ( 1995) “Oscillating flow about two and three dimensional bilge keels”, Proc. 14th Int'l. Conf. of Offshore Mechanics & Arctic Engineering (OMAE), Copenhagen, Denmark.
[19] Schmitke, R.T. ( 1978) “Ship sway, roll, and yaw motions in oblique seas”, Trans. SNAME, 86, pp. 26–46.
[20] Singh, S. ( 1979). “Forces on bodies in oscillatory flow”, Ph.D. dissertation, Univ. of London, UK.
[21] Susbielles, G. and Bratu, C. ( 1981). vagues et ouvrages pétroliers en mer, Publications of The Institute of Petroleum of France, Editions Technip, Paris, France.
[22] Tanaka, N. ( 1960). “A study on the bilge keels. Part 4. On the eddymaking resistance to the rolling of a ship hull”, J. Soc. Naval Arch. Japan, 109.
[23] Vaidhyanathan, M. ( 1993)., “Separated flows near a free surface”, Ph.D. dissertation, Dept. of Naval Arch. and Offshore Engrg., Univ. of Calif., Berkeley.
[24] Vinje, T. and Brevig, P. ( 1980). “Nonlinear, twodimensional ship motions” Tech. Rept., The Norwegian Inst. of Tech., Trondheim.
[25] Vugts, J.H. ( 1968). “The hydrodynamic coefficients for swaying, heaving and rolling cylinders in a free surface”, Rept. No. 194, Laboratorium voor Scheepsboukunde, Technische Hogeschool Delft, The Netherlands.
[26] Wehausen, J.V. and Laitone, E.V. ( 1960). “Surface Waves”, in Handbuch der Physik, 9, SpringerVerlag.
[27] Willert, C.E. and Gharib, M. ( 1991). “Digital particle image velocimetry”, Experiments in Fluids, 10, pp. 181–193.
[28] Yeung, R.W. ( 1982). “Numerical methods for freesurface flows”, Ann. Rev. Fluid Mech., 14, pp. 395–442.
[29] Yeung, R.W. and Ananthakrishnan, P. ( 1992a). “Oscillation of a floating body in a viscous fluid”, J. Engrg. Math., 26, pp. 211–230.
[30] Yeung, R.W. and Ananthakrishnan, P. ( 1992b). “Vortical flows with and without a surfacepiercing body”, Proc. 19th Symp. on Naval Hydrodyn., Seoul, Korea.
[31] Yeung, R.W., and Cermelli, C.A. ( 1996), “Vortical Flow Generated by A Plate Rolling in a Free Surface”, Chapter in Advances in Fluid Mechanics, Computational Mechanics Publisher.
[32] Yeung, R.W., Sphaier, S.H., and Vaidhyanathan, M. ( 1993). “Unsteady flow about bluff cylinders”, Int'l. J. Offshore and Polar Engrg., 3, pp. 81–92.
[33] Yeung, R.W. and Vaidhyanathan, M. ( 1993). “Flow past oscillating cylinders”, ASME J. Offshore Mech. and Arctic Engrg., 115, pp. 197–205.
[34] Yeung, R.W. and Vaidhyanathan, M. ( 1994). “Highly separated flows near a free surface”, Proc. Int'l. Conference on Hydrodynamics, Wuxi, China.