The prediction of roll motion of a ship with bilge keels is particularly difficult because of the nonlinear characteristics of the viscous roll damping. Flow separation and vortex shedding caused by bilge keels significantly affect the roll damping and hence the magnitude of the roll response. To predict the ship motion, the Slender-Ship Free-Surface Random-Vortex Method (SSFSRVM) was employed. It is a fast discrete-vortex free-surface viscous-flow solver developed to run on a standard desktop computer. It features a quasi-three-dimensional formulation that allows the decomposition of the three-dimensional ship-hull problem into a series of two-dimensional computational planes, in which the two-dimensional free-surface Navier–Stokes solver Free-Surface Random-Vortex Method (FSRVM) can be applied. In this paper, the effectiveness of SSFSRVM modeling is examined by comparing the time histories of free roll-decay motion resulting from simulations and from experimental measurements. Furthermore, the detailed two-dimensional vorticity distribution near a bilge keel obtained from the numerical model will also be compared with the existing experimental Digital Particle Image Velocimetry (DPIV) images. Next, we will report, based on the time-domain simulation of the coupled hull and fluid motion, how the roll-decay coefficients and the flow field are altered by the span of the bilge keels. Plots of vorticity contour and vorticity isosurface along the three-dimensional hull will be presented to reveal the motion of fluid particles and vortex filaments near the keels.

## Introduction

Consideration of the static and dynamic stability properties of a ship is fundamental to its design and safe operation. The stability of a ship, which is needed to avoid a ship from capsizing, is directly connected to its roll motion. The natural period of roll is designed to be away from the period of the higher energy of the wave spectrum. However, the encountering frequency and direction of waves could lead to the coalescing of the two periods, which may cause capsizing or severe damages because of the resulting large motion. In the unfortunate situation that resonant roll motion takes place, large damping is the only recourse to reduce the hazardous response. Bilge keels, as shown in Fig. 1, have been the traditional passive “stability enhancement system,” offering an increase in the hydrodynamic resistance when a ship rolls, thus limiting the amount of roll a ship has to endure and yet commanding little increase in operation costs. The primary damping mechanism arising from a bilge keel is the formation and shedding of vortices. These vortices significantly affect the roll damping and the evolution of the roll response.

To predict the roll motion of a hull, a number of theoretical and experimental approaches have been taken (e.g., see Refs. [1–4]). Inviscid-fluid numerical methods (including inviscid strip theories, panel methods, and system-based filtering methods) have been widely used in many different applications, but they usually require empirical or experimental data to estimate the linear and nonlinear roll damping coefficients. For various types of applications, the empirical nonlinear damping may rely on different formulations, with the aim to solve the ordinary differential equation for roll (e.g., see Refs. [5,6]). It is known that different types of nonlinearities in a nonlinear system can lead to diverse characteristics in terms of the system response. Identifying the right type of nonlinear roll damping is also a difficult task. However, the dynamic interaction of fluid flow and roll motion may also be obtainable from first principles. Recently, large-scale computational fluid dynamic (CFD) method has also been applied in this area and shows promise (e.g., see Refs. [7,8]), but it may require finer elements to achieve high-fidelity simulation. Issues of mesh density and long computational time often prevail, and the CFD method may not be the most effective analysis tool in a design process with multiple parameters.

In parallel to classical inviscid theories and CFD, the Free-Surface Random-Vortex Method (FSRVM) has also been developed and applied to simulate flows near sharp keels and fins (e.g., see Refs. [9–12]). These earlier developments, however, were aimed at simulating two-dimensional floating bodies.

To predict the roll motion of a three-dimensional vessel with bilge keels, we have completed some recent theoretical development, using the UC-Berkeley code, called Slender-Ship Free-Surface Random-Vortex Method (SSFSRVM). It is a fast free-surface hydrodynamics solver designed to run on a standard desktop computer and features a quasi-three-dimensional formulation involving decomposing the problem into a series of two-dimensional computational planes. The theory for this three-dimensional to two-dimensional conversion was conceived in Refs. [13] and [14] for inviscid fluid. However, its extension to a viscous fluid was not straightforward and was obtained using scaled variable analysis by Yeung et al. [15], which allows the two-dimensional computational engine of FSRVM [16] to be used for three-dimensional flow. Previously, the most advanced applications were given in Ref. [17], in which the exact body boundary condition was satisfied on the instantaneous wetted surface of the moving vessel with inviscid but fully nonlinear free-surface boundary conditions. Documentation for these forced or prescribed motion cases is given in Refs. [15,17]. SSFSRVM has recently been developed to accommodate the free-body motion in the time domain to study the roll motion of ship hulls. This is carried out in a manner similar to the two-dimensional treatment of Roddier et al. [18], allowing one to model the dynamic coupling between body motion and fluid motion, but now with the full capability of four degrees-of-freedom (DOF) motion: free sway, heave, roll, and pitch motions in three dimensions. The latest version of SSFSRVM model does not require the assumption of small amplitude motion and is capable of having viscosity turned on or off in the solution procedure. Because FSRVM uses a mesh-free formulation, there is no issue with numerical viscosity, and the method is efficient in producing accurate predictions at a fraction of the time required by methods such as Reynolds-averaged Navier-Stokes (RANS) [19].

On the other hand, SSFSRVM, similar to other methods, has some limitations. For instance, this method can be only applied to simulate slender bodies. Unlike the mesh-based CFD methods, it only computes the longitudinal vorticity component; therefore, vortex stretching and tilting are not accounted for [15]. Also, the attached turbulent boundary layer is not fully resolved. This method is designed to capture the dominant vortical structures generated by the bilge keels and predict its effect on the ship motion as efficiently as possible.

In this paper, the effectiveness of SSFSRVM modeling will be examined first by a comparison of the time histories of roll decay resulting from simulations and experimental measurements. Furthermore, the evolution of vorticity contour near the port-side bilge keel will be presented to validate the accuracy of the model. We will also report how the free roll motion and the roll decay coefficients are quantitatively altered by the span of the bilge keels. Finally, the roll responses of the hulls with different bilge keels in regular incident beam waves will be presented and discussed.

## Theoretical Model

Based on the slender body theory, a three-dimensional problem can be reduced into a sequence of two-dimensional problems, for an adequately slender hull. The approximation is not obvious if the fluid is viscous and is discussed in Ref. [15]. The hull is geometrically defined by a series of sectional profiles (or stations) which are equally spaced along the length of the vessel, see Fig. 2.

For the case with zero forward speed, the slender body theory is essentially the same as the strip theory. Two-dimensional computational planes are set up along the hull and fixed at their longitudinal positions on earth to simulate the disturbance generated by the unsteady motion of the three-dimensional hull. For the case with a forward speed, the two-dimensional computational planes are still fixed on earth. However, because of the forward motion of the ship, the two-dimensional body contour cut by each computational plane is deforming in time. The concepts of pseudo time and expansion velocity are adopted to modify the body boundary condition (see Ref. [15]). In addition, at each time step, a new two-dimensional computational plane is added at the first bow section, and the last computational plane passed by the stern is removed from the computational procedure.

The transient flow problem in each two-dimensional cross-plane is solved by FSRVM, which is an efficient two-dimensional free-surface flow solver coupled with the random vortex method of Chorin [20] to account for viscous effects. It relies on a vorticity and stream-function formulation, coupled with nonlinear free-surface and body boundary conditions for modeling wave–body interactions [21]. The solution is obtained by decomposing the flow field into an irrotational component and a vortical component. The irrotational component of the flow is solved using a complex-variable Cauchy integral method, based on the instantaneous geometry of the computational domain and the vorticity field at that instant. The vorticity field is solved using the random-vortex method for the diffusion effects and an order N multipole-accelerator interaction algorithm for convection effects. More details can be found in Ref. [10].

The theoretical part of this paper will focus on the derivation of the 4DOF free-body equations of motion and their solution by taking advantage of FSRVM computational engine.

### Coordinate Systems.

For a ship hull, six independent coordinates are necessary to represent its position and attitude: surge, sway, heave, roll, pitch, and yaw. In the current SSFSRVM model, three different coordinate systems are used to develop and solve the rigid-body motion of the vehicle (see Fig. 3). The first one is an earth-fixed coordinate system $O\xafx\xafy\xaf\chi \xaf$, which is used to record the translational motion and the rotation of the vessel. The $(O\xafx\xaf\chi \xaf)$-plane lies on the still water surface, the positive $\chi \xaf$-axis is in the direction of the forward speed, and the positive $y\xaf$-axis points upwards. We use *χ* instead of the conventional *x* to denote the forward direction as it is desirable to name the sectional-plane variables $(x,y)$ so that the complex variable $z=x+iy$ can be used as in FSRVM [10]. At the initial time, the origin of the earth-fixed coordinate system coincides with the geometrical center of the ship, $O\u0302$, which is the intersection point of the still water plane, the midship plane, and the centerplane.

The second coordinate system is a body-fixed coordinate system $O\u0302x\u0302y\u0302\chi \u0302$, which is used to define the position of the body nodes, with $O\u0302$ attached to the geometrical center of the ship. $\chi \u0302$ points in the forward longitudinal direction of the ship, $x\u0302$ in the port-side direction, and $y\u0302$ in the upward direction.

The third coordinate system $Oxy\chi $ is a steadily moving frame of reference that has the same constant forward speed as the vessel. When the vessel does not have a forward speed, $Oxy\chi $ would be fixed on earth and is essentially the same as the earth-fixed coordinate system $O\xafx\xafy\xaf\chi \xaf$. There are N translating coordinate systems with one for each two-dimensional computational plane. Each one of these translating sub-systems is set up to solve the two-dimensional transient flow problem. Figure 4 illustrates an example of the translating *sub-systems*, with origin *O _{i}*. The directions of the $xiyi\chi i$ axes of each sub-system are set to be the same as those of the translating $xy\chi $ axes. Therefore, the $(Oixi\chi i)$-plane also lies in the still water surface with $yi$-axis pointing upward. Among the translating coordinate systems, the one at mid-ship is chosen to describe the global equations of motion of the vessel and is named as

**referenced**translating coordinate system, $Oxy\chi $. The others are named as translating sub-coordinate systems. At the initial time, the origin $O$ is located at the origin of the earth-fixed coordinate system. Since the vessel could be allowed to move forward with sway and heave motion, the body center, $O\u0302$, would not necessarily be coincident with the referenced origin

*O*.

### Relation Between Global Body Motion and Sectional Motion.

Figure 5 shows an example of the FSRVM model in a two-dimensional computational plane. The *i*th section of a vessel translates and rotates in this fluid domain according to the global motion of the point $O\u0302$. Point $Oi\u0302$ is the geometrical center of the local section. Its location with respect to the *i*th translating subframe is denoted by $(xO\u0302i,yO\u0302i)$.

*i*th section, with initial coordinates $(x\u0302bo,i,y\u0302bo,i)$ on the body, will have new coordinates $(xbo,i,ybo,i)$ in the $Oixiyi$ system:

where $(xO\u0302i,yO\u0302i)$ can be obtained from the motion of the reference point $O\u0302$ of the vessel, and *θ _{b}* is the Euler roll angle about the

*χ*-axis.

The hull section in each computational plane has, evidently, only 3DOF (sway, heave and roll). Since the actual vessel has 6DOF, the sectional motion of sway, heave and roll is directly related to the global motion of the $\u0302$ system. As the influence of surge motion (in the presence of forward velocity) is modeled by the concept of expansion velocity (see Ref. [15]), global pitch and yaw motion leads to sectional sway and heave motion. The following can be deduced in a straightforward manner: per Fig. 6, a positive (counter-clockwise from above) yaw of the vessel, $\alpha b$, induces sway of sections from bow to stern, while a positive (clockwise from starboard side) global pitch, $\gamma b$, induces heave of sections from bow to stern. Hence, with the contribution of global yaw and pitch, the position and rotation of the *i*th sectional center $(xO\u0302i,yO\u0302i)$ can be expressed in the following form:

where $xb,\u2009yb,\u2009\theta b,\u2009\gamma b$, and $\alpha b$ denote the sway, heave, roll, pitch, and yaw displacements of the vessel with respect to the translating frame. $Li$ denotes the distance from the *i*th section to the midship section, with bow sections positive and stern sections negative.

### Equations of Motion.

In order to simulate the free response of the floating body, the equations of motion must be written and solved with respect to one reference point in a coordinate system. Normally, equations of motion are written about the center of gravity (COG) of a body. Alternatively, it is preferable to develop these equations at the body center $O\u0302$, rather than *G*. In this section, it is more convenient to obtain the equations of motion with respect to point *G* first (see Fig. 3) and then transform the results to $O\u0302$ in the steadily translating coordinate system.

where *M _{b}* is the body mass, ($\chi \u0308g,\u2009x\u0308g,\u2009y\u0308g$), and ($\theta \u0308g,\u2009\gamma \u0308g,\u2009\alpha \u0308g$) denotes the translational and angular accelerations (the Euler's angles) of the body. The tilde $\u2009\u0303$ over each variable in Eq. (4) denotes that the value of this variable is calculated about point

*G*in the translating coordinate system.

*G*can be obtained from the moments measured at point $O\u0302$ by

where $FG$ and $MG$ denote the force and moment vectors about point *G*, respectively; $MO\u0302$ denotes the moment vector about point $O\u0302$; and $RG$ denotes the vector from $O\u0302$ to *G* with respect to the translating frame. $RG\u2261(\chi \u02c7g,x\u02c7g,y\u02c7g)T=(\chi g\u2212\chi b,xg\u2212xb,yg\u2212yb)T$. The $\u02c7$ over each variable denotes that the value of this variable is calculated about point $O\u0302$ with respect to the translating coordinate system.

*G*and $O\u0302$ will have the same values. So,

*G*, $aG$, can be expressed in terms of the acceleration at point $O\u0302$, $ab$, by the following equation:

where $I\chi \chi $ and $Ixx$ denote the roll and pitch moments of inertia about point $O\u0302$, respectively.

### Total and Sectional Forces and Moments.

*i*th station, and $F3,i$ denotes the sectional roll moment. Apart from the first equation, showing the limits of the summation over

*i*, we will henceforth adopt the simplified notation of $\u2211i$ to mean the same limits of

*i*. The sectional loads can be obtained by integrating the pressure over the wetted sectional body contour $\u2202Db$,

where *p* denotes the fluid pressure that can be calculated by using Euler's integral.

*i*th section with respect to the corresponding translating frame of reference as follows:

where $Ajk,i,\u2009j,k=1,2,3$, represent the sectional hydrodynamic added mass or added moment of inertia; $A4k,i,\u2009k=1,2,3$, represent the summation of sectional hydrodynamic damping and hydrostatic restoring loads; and $\Delta L$ denotes the distance between any two adjacent sections.

*B*and

_{jk}*W*are defined as:

_{jk}Equation (16) completely describes the full nonlinear dynamic coupling between the fluid and three-dimensional vessel, which can be solved explicitly to yield $(x\u0308b,y\u0308b,\theta \u0308b,\gamma \u0308b)$ at any given time, *t*. The coupling with the FSRVM code occurs in the fluid pressure *p* computations, and the accelerations solved determine the body boundary condition of the fluid problem.

Further, with the acceleration being solved, the velocity and position can be integrated in time to advance the location of the body. Then, the kinematic and dynamic conditions are used to update the boundary nodes. A new set of boundary conditions is available for the next time step. The time-domain simulation for free motion can now be achieved.

## Effectiveness of Slender-Ship Free-Surface Random-Vortex Method Modeling

Aloisio and Di Felice [22] conducted an experimental study of the velocity field around a ship model in free roll decay motion at the INSEAN towing tank No. 2. The flow around the bilge keels was investigated using a two-dimensional particle image velocimetry (PIV) underwater system. Their experiments visualize the fluid dynamics and provide a set of certified experimental data for validating our numerical method. To evaluate the accuracy of SSFSRVM in predicting the roll motion of a three-dimensional hull with bilge keels, these experimental results are compared with our numerical results.

In the experiments, a 5720 mm length model (INSEAN C2340) with bilge keels (see Fig. 1) is chosen and tested in the free roll decay motion. Measurements were performed with an initial roll angle of 10 deg. The span of the bilge keel is 4.76% of the full beam of the hull.

### Free Roll Decay Response at Zero Forward Speed.

By using SSFSRVM viscous-flow model, we simulated the same roll-decay test at zero forward speed carried out by Aloisio and Di Felice [22]. The comparison of the time histories of roll motion from the experiments and SSFSRVM simulations is presented in Fig. 7. The numerical solution contains full memory effects from the free surface and, of course, from the vorticity field. It can be seen that the predicted time history matches very well with the experimental data, a remarkable confirmation of the present theory that is based on first principles.

To further validate the accuracy of the model on simulating the viscous-fluid flow around the ship for the zero speed case, comparisons of the vorticity evolution near the bilge keel on the port-side (at the longitudinal position, $\chi /L=0.504$) between PIV measurements and simulations are presented in Figs. 8 and 9. The blue color denotes negative (counter-clockwise, when observed from stern) vorticity, while the red color denotes positive (clockwise) vorticity. The predicted vorticity field uses the same color-bar scale as the PIV measurements. It is observed, from both PIV measurements and simulation, a vortex rolls up gradually at the keel with increasing strength and core size during the first half period of oscillation. A pair of counter-rotating vortices is generated as the roll motion changes directions. The distances that the vortex pair travels are also predicted closely by the numerical method, as shown in Fig. 9. Since SSFSRVM only simulates the longitudinal vorticity field, minor differences can be found between the simulated and the experimental results. Overall, we consider the model to be capturing the dominant motion of the vortical structures well.

The vorticity evolution in Figs. 8 and 9 illustrates the motion of vortices and validates the numerical model in a sectional plane. To visualize the predicted three-dimensional nature of the vortical structures, Fig. 10 presents a time sequence of vorticity contours and vorticity iso-surfaces along the three-dimensional hull in the first cycle of roll decay motion. Times are mainly selected to correspond with the peaks and the zero crossing of the roll angle as shown in the upper left legends of Fig. 10. In these figures, the blue and red colors still denote the negative and positive vorticities, respectively. In the lower right corner, a two-dimensional vorticity contour at the midship section is displayed as a sectional diagram of the three-dimensional structures. The iso-surfaces are plotted in the three-dimensional sub-figures to represent the vortex filaments in the fluid. In the following discussion, we will focus on the starboard bilge keel for convenience.

Figure 10(a) shows the vorticity distribution at the moment when the hull starts to change its roll direction. It can be seen that a positive strong vortex filament has rolled up with strong vortex cores along the starboard bilge keel. Besides, a negative vortex filament is just generated at the tip of the keel with a very small core. As the roll motion progresses in time, the newly generated negative vortex filament rolls up. When the ship is in the upright position (see Fig. 10(b)), a fully formed vortex-filament pair is evident. Then, the interaction between the counter-rotating vortices convects the pair away from the hull and rapidly mixes their energy, as shown in Fig. 10(c). When the ship reverses its direction (see Fig. 10(d)), a new negative vortex filament is generated near the keel tip. At this moment, the previously generated vortex pair has split into parts and almost dissipated into the surrounding fluid. As the roll motion progresses in time, a positive vortex filament will be generated at the keel. Then, a new cycle of vortex shedding will start.

These four figures illustrate the motion of vortex filaments in the first oscillation. We find that bilge keels alternately generate positive and negative vortex filaments and form vortex pairs along the hull. The vortices, which represent high energy dissipation from the hull's roll motion, increase the viscous damping and reduce the roll amplitude. These figures reveal the working mechanism of the bilge keel from the energy-dissipation point of view. Through the earlier analysis, we consider that the SSFSRVM model is capable of simulating the dominant three-dimensional vortical structures.

In terms of the computation time, SSFSRVM takes five clock hours on a standard desktop computer with a four-core central processing unit (CPU) to simulate ten cycles of roll motion of a ship model (5.7 m long) at zero forward speed. The CPU time per cycle of roll motion is about 2 h.

Kristiansen et al. conducted a numerical study on the effect of bilge keels on a ship model with a similar length by using an RANS method [23]. They also numerically simulated the free decay with the same initial angle of 10 deg. Regarding the computational time of the RANS method, running the eight periods took 34 h on 128 cores in parallel. The corresponding CPU time is 544 h per cycle of roll motion. The difference in CPU time clearly reveals the computational efficiency of the SSFSRVM.

Figure 11 presents the free decay test of the vessel in 4DOFs with an initial roll angle of $15deg$. The vessel was given zero initial heave, pitch, and sway motions. The coupling between different motion modes transfers energy from roll motion to heave, pitch, and sway as well, resulting in an intricate 4DOF motion.

### Free Roll Decay Response With Forward Speed.

In addition to the free decay test at zero forward speed, Aloisio and Di Felice [22] also conducted a decay test at the same INSEAN towing tank with a model scale forward speed of 1.03 m/s, corresponding to the Froude number $Fr=0.138$. In this section, we examine the accuracy of SSFSRVM model in simulating the roll motion of the ship model at forward speed.

In the experiment, the ship model was initially towed in calm water with the target forward speed in static heel to port. When the surface waves were in a steady-state, the ship model was released from an initial angle of 10 deg by using an electromagnet switch. The same sequence of events is simulated by SSFSRVM. The numerical ship model pierces in calm water with a constant hull speed of 1.03 m/s and with a constant roll angle of $\theta 0=10deg$ for about 5 s when the pattern of the surface waves does not change. Then, the ship is released and allowed to roll freely from the $10deg$. The comparison of the time histories of roll motion for this scenario resulting from the experiments and SSFSRVM simulations is presented in Fig. 12. Again, the predicted time history closely matches the experimental data, especially for the first three periods. In addition, the predicted bilge keel vorticity evolution at the longitudinal location $\chi =\u22120.175L$ has been compared with the experimental DPIV images, as shown in Fig. 13. The same level of agreement has been found as the zero speed case.

Based on the previous comparisons of roll time histories and vorticity evolution, we validate the SSFSRVM in simulating the roll motion of a vessel with and without forward speed.

## Effects of Bilge Keels

Four SSFSRVM models with different configurations of bilge keels are built to evaluate the effects of the bilge keels, as shown in Fig. 14. These four models have the same hull design (INSEAN model C2340) but different bilge keel spans. The span of BK1 is 4.76% of the full beam of the hull, which is the same as that of the model C2340. The spans of BK2 and BK3 are two and three times of the span of BK1, respectively, while the model BK0 has no bilge keel. Figure 15 shows the midship sections of these four models to provide more details.

### Roll Decay Coefficient.

The free roll decay tests of these four models were conducted with the same initial roll angle of 10 deg. The time histories of roll motion and the hydrodynamic normal forces acting on the port-side bilge keel are shown in Figs. 16 and 17. From these figures, it can be seen that a larger bilge keel provides larger damping to decrease roll amplitude and also provides larger hydrodynamic added moment of inertia to increase the roll period. At the same time, the larger bilge keel experiences much larger hydrodynamic loads. In Fig. 17, increasing the bilge keel span dramatically amplifies the normal force. Though the large keel reduces the amplitude of motion of the ship, the increased normal force may threaten the safety of the hull structure.

*B*

_{1}and quadratic damping

*B*

_{2}is given by the following normalized equation:

*ω*is the natural frequency, and

_{n}*B*

_{1}and

*B*

_{2}are often obtained by experimental means or regression analysis. The SSFSRVM procedure can provide these coefficients from solution of the roll response. If energy principle is applied to this simple oscillator, it can be shown that $n$, defined as the logarithmic decrement, is a function of mean roll amplitude $\eta \xaf$ (see Ref. [24])

where *η _{i}* and $\eta i+1$ correspond to successive double roll-amplitudes (see Fig. 18). These two measurable quantities can be used to deduce the damping coefficients:

Thus, *B*_{2} and *B*_{1} will show up as the slope and the *y*-axis intercept of $n$, respectively, when $n$ is plotted versus $\eta \xaf$. Under these assumptions, the nonlinear characteristics of roll damping can be characterized by such $n$ versus $\eta \xaf$ plots. Deviations from the relation of Eq. (20) implies the imperfect modeling of the complex fluid–structure interaction by such a simple two-term representation.

Figure 19 displays the non-dimensional decay coefficient versus the mean roll angle derived from the roll decay data (Fig. 16). In the experimental roll decay test, the ship model is released from an initial roll angle in a fluid at rest; therefore, there is no existing flow structure prior to the start of the roll motion. Our numerical investigation follows the same procedure. The memory effect will alter the decay coefficients obtained from the first few roll periods (see Ref. [25]). Therefore, when we calculated the linear regression of the decay coefficient, the first three roll periods with roll amplitude larger than 6 deg are excluded from the damping coefficient calculations. The linear trend lines of the valid decay coefficients are shown in Fig. 19 for all the four models.

It is observed that the trend-lines fit the computed data points fairly well. This indicates that the quadratic component of the total damping plays an important role when the mean roll angle is small. As expected, the BK3 model offers the largest quadratic damping coefficient, since its trend-line has the largest slope. In addition, a monotonic increase of the decay coefficient with the bilge keel span can be seen in the figure. However, the increase of the decay coefficient from the BK1 to the BK2 model is larger than that from BK2 to BK3. Thus, the quadratic roll damping does not increase linearly with the span of the bilge keel.

To examine the differences in the roll decay coefficient *n* of these four models when the roll angle is small, the vorticity contours at the midship section in the seventh oscillation obtained from SSFSRVM simulations are presented in Fig. 20. In these figures, each column represents one bilge keel size. By comparing plots in each row, we can examine the effects of bilge keels.

In this cycle of roll, the motion amplitude is about $3deg$. From this figure, we note that for the model without bilge keels, the vorticity stays close to the body. That is, there is almost no flow separation, and the surrounding fluid particles have little vortical motion. However, for the other models, it can be seen that strong vortical structures are generated by the bilge keels, particularly the hulls of BK2 and BK3. It is known that vortical structures characterize high energy dissipation. By comparing the vorticity contours of BK2 and BK3, we find that the vortices generated by these two models do not have significant difference with respect to their strength and core size. However, the vorticity field of the BK3 model appears to be more complex and dynamic. This indicates that slightly more energy is dissipated into the fluids by the larger bilge keel and explains the fact that BK3 model has a little larger decay coefficient than the BK2 model in Fig. 19.

### Roll Response in Beam Waves.

Slender-ship free-surface random-vortex method can be also applied to simulate the roll response of a ship in beam waves. To demonstrate this, an oscillating pressure patch is used to generate regular incident waves on the starboard-side free surface with the wave height chosen to be 2.48 cm, which is $10%$ of the draft of the ship hull. The wave period is 2.25 s. In beam waves, the ship can only roll freely in waves with the other motions constrained.

Figure 21 shows the time histories of the roll responses among hulls with different bilge keel spans. Since the wave period is close to the roll resonance period of these vessels, the roll motion ramps up rapidly, and the motion amplitudes are quite large in the steady-state. Through the comparison of the responses among these four models, we find that the roll amplitude in the steady-state decreases with the increase of the bilge keel size. The largest bilge keels reduce the motion amplitude most significantly. Compared to the BK0 model, the motion amplitudes of the BK1, BK2, and BK3 models are reduced by 6%, 11%, and 32%, respectively.

## Conclusions

In this study, we proposed and have successfully implemented the SSFSRVM in simulating the free roll decay motion of a hull. The predicted time histories of the roll motion and the vorticity distribution near a bilge keel are compared with the experimental measurements and the agreements are very good. The validation works demonstrated that the SSFSRVM model is capable of predicting the roll motion of a hull and simulating the dominant behavior of the longitudinal vortical structures in the fluid, in comparison with documented PIV images. With respect to the computational efficiency, the SSFSRVM model is much more efficient than mesh-based CFD methods and yet can capture significant flow details.

With regard to the effects of bilge keel span, it is found that a larger bilge keel generates stronger counter-rotating vorticity pairs. This transfers more energy from the hull into the surrounding flow and leads to a larger decay rate. The roll decay coefficient increases with the bilge keel span, but this rate of increase is not proportional to the rate of increase of the bilge keel span. In the cases of roll motion in waves, we found that increasing the bilge keel span effectively decreases the roll amplitude.

Overall, considering the accuracy and the computational efficiency, we believe that SSFSRVM is an excellent modeling and simulation tool for predicting the motion of a ship with bilge keels. This method provides an effective tool to evaluate and improve the hull performance in modern ship design.

## Acknowledgment

The authors would like to acknowledge Dr. G. Aloisio and Dr. F. Di Felice of INSEAN for their high-quality experimental work. We are grateful to Lu Wang of the Berkeley Marine Mechanics Laboratory for his valuable assistance. Partial support for the work in this paper was enabled by the American Bureau of Shipping Endowed Chair in Ocean Engineering (RWY).

### Appendix Convergence of Numerical Method

The numerical errors and uncertainties caused by the boundary element size and time step are examined. The INSEAN ship model C2340 with bilge keels is chosen and built numerically. The three-dimensional ship model is decomposed into 51 two-dimensional planes.

#### Boundary Element Size

Three different element sizes are used to generate the boundary of the computational domain in order to estimate numerical errors. The coarse model has 200 elements on each two-dimensional plane, the medium-element model has 282, and the fine model has 400. The corresponding element size normalized by the full beam of the hull $Sg/B$ is around 0.0280, 0.0198, and 0.0140 for the coarse, medium, and fine elements, respectively. A refinement ratio of $2$ is used in this study. A small time step $\u25b3t=0.0115$ s is selected to make sure that the simulations are convergent.

Figure 22 shows the comparison of the time history of the roll-decay motion at zero forward speed among different models. It can be seen that the element size does not affect the result much, especially during the first several oscillations. Figure 23 shows the solution difference between the coarse/medium elements, $\epsilon g12=\theta g1\u2212\theta g2$, and the medium/fine elements, $\epsilon g23=\theta g2\u2212\theta g3$. Examination of $\epsilon g23$ and $\epsilon g12$ shows a monotonic convergence $\epsilon g23<\epsilon g12$. The maximum difference of the roll response for the medium versus fine elements is less than 0.12 deg. Considering the accuracy and efficiency, we choose the medium element in the present study.

#### Time Step Size

The same procedure is conducted to estimate the temporal discretization error. The refinement ratio is $rT=\u25b3t2/\u25b3t1=\u25b3t3/\u25b3t2=2$. The coarse, medium, and fine time steps are selected as $\u25b3t1=0.046$ s, $\u25b3t2=0.023$ s, and $\u25b3t3=0.0115$ s, respectively. Correspondingly, each period of the roll decay motion is roughly resolved with 50, 100, and 200 time steps. The medium element is used for all the three cases.

Figure 24 shows a comparison of the time history of the roll motion among different time step models. It can be seen that the solution has converged. The difference between the medium/fine time steps is almost not noticeable. Again, considering both the computational accuracy and efficiency, we choose the medium time step as the time step in this work.

In summary, we find that the solution of SSFSRVM is not sensitive to the element size and time step. The uncertainty level is low. The combination of 300 elements on each section and 100 time steps per roll period *T* is a good choice considering the balance between accuracy and computational efficiency.