Near field flow characteristics around catamarans close to resonant conditions involve violent viscous flow such as energetic vortex shedding and steep wave making. This paper presents a systematic and comprehensive numerical investigation of these phenomena at various oscillating frequencies and separation distances of twin sections. The numerical model is based on the solution of Navier–Stokes equations assuming laminar-flow conditions with a volume of fluid (VOF) approach which has proven to be particularly effective in predicting strongly nonlinear radiated waves which directly affect the magnitude of the hydrodynamic forces around resonant frequencies. Considered nonlinear effects include wave breaking, vortex shedding and wave-body wave-wave interactions. The method is first validated using available experiments on twin circular sections: the agreement in a very wide frequency range is improved over traditional linear potential flow based solutions. Particular attention is given to the prediction of added mass and damping coefficients at resonant conditions where linear potential flow methods fail, if empirical viscous corrections are not included. The results of the systematic investigation show for the first time how the so-called piston-mode motion characteristics are nonlinearly dependent on the gap width and motion amplitude. At low oscillation amplitudes, flow velocity reduces and so does the energy lost for viscous effects. On the other hand for higher oscillation amplitude, the internal free surface breaks dissipating energy hence reducing the piston mode amplitude. These effects cannot be numerically demonstrated without a computational technique able to capture free surface nonlinearity and viscous effects.

## Introduction

The discussion on the resonant phenomena in oscillating twin sections often goes back to the work of Molin [1] who investigated the so-called piston mode resonant fluid motion inside moonpools in large rectangular barges. The piston mode motion of the fluid inside twin hulls consists in the stationary oscillation of a water column whose amplitude is enhanced with respect to the outer radiated wave. Resonant phenomena have been deeply investigated by McIver et al. [2], who proposed analytic and numerical solutions for the trapping mode and resonant mode motions between surface-piercing elements in two dimensions (see Refs. [3,4]) and in three dimensions using a high order linear panel method (see Ref. [5]).

Fluid dynamic research has been focusing on this topic for many years: its relevance can be in fact generally related to trapped water between two vertical walls, e.g., moonpools or ship moored by a bottom mounted terminal (see Ref. [6]).

Many studies have been presented proposing different techniques for the prediction of resonance frequencies and water motion amplitude occurring at these particular oscillation modes. One of the first studies on the interaction between the inner and the outer free surfaces was presented by Molin using a potential flow theory [7]. A linear free surface condition has been included in a method solving for the velocity potential in the domain laterally bounded by the two inner walls and vertically bounded by the bottoms of the two hulls. Molin modeled water flow around the keels using two sinks, some distance away from the hull. This method neglects the existence of radiated waves, hence causing a phase shift of the fluid motion inside the moonpool with respect to the forced body motion. The analytic solution of radiated waves problem has been addressed using a domain decomposition technique by Faltinsen et al. [8], Yeung and Seah [9] and Porter and Evans [10] with methods still relying on potential flow assumption. These techniques are yet based on the small amplitude hypothesis and potential flow assumption with a linear free surface condition. Nonlinear and viscous effects are known to significantly affect the radiation problem for twin oscillating bodies in calm water (see Ref. [11]). Linear potential flow based theories demonstrate weaknesses when applied to resonant gap problems. Viscous and nonlinear effects play in fact an important role in damping out the very large water motions occurring at resonance: potential flow theories may result in drastic failure if not properly corrected. To address this problem, Timokha and Faltinsen [12] quantified theoretically the pressure discharge in a two-dimensional (2D) moonpool opening based on a similarity between piston-mode sloshing in a two-dimensional rectangular moonpool and a small-amplitude oscillatory flow through a slot in a slatted screen. Buchner et al. [13] implemented a lid boundary condition at the gap free surface with the aim to suppress the nonrealistic high wave elevation resulting from linear potential flow models. Also in this case, a linear damping coefficient resulting from experimental decay tests is required. The lid method was also used by Newman [14] and Chen et al. [15] who reformulated the traditional free surface condition using an empirical parameter *ϵ* found from computational fluid dynamics (CFD) or experiments. The improvements brought in by the lid method have been demonstrated also by Pauw et al. [16] who concluded that an empirical tuning of damping coefficient is necessary since an unique value cannot be found. Centeno et al. [17] proposed a numerical method based on a potential flow approach where viscous effects were recovered through Thwaites cross flow theory; however, this simplified theory cannot reach anything more than a qualitative agreement (see Ref. [18]).

Kristiansen and Faltinsen [6] confirmed the limits of linear potential flow theories applied to gap resonance problem: especially for piston-mode motions the free surface becomes extremely nonlinear and the assumption of small amplitude motion can no longer be applied. Viscous damping and vortex shedding can explain the discrepancies near resonance.

Kristiansen and Faltinsen applied a boundary element method with mixed Eulerian–Lagragian approach for the interface tracking where the effects of flow separation have been recovered using an inviscid vortex tracking method. They concluded that including vortex shedding improves numerical results, demonstrating the weakness of traditional boundary element method (not including separation) when applied to this problem.

Lu et al. [19] studied the effect of incident waves on a multihull configuration (three fixed bodies) using both a linear potential method and a finite elements viscous solver based on a volume of fluid (VOF) technique for the free surface capturing. Recently, Kistiansen and Faltinsen [20] have proposed a very effective domain-decomposition method combining potential and viscous flow. The potential flow solver used to propagate waves in the far field allows for an improved computational efficiency and, if compared with surface capturing technique, it also allows for a more accurate solution of small amplitude waves radiated at a large distance. The domain of the viscous flow solver is limited to the near field region. This method has been recently validated in three-dimensional by Kristiansen et al. [21], successively adapted to study forward speed effects by Fredriksen et al. [22], and adopted by Fredriksen et al. [23] for the solution of incident waves on a floating two-dimensional body.

It is straightforward to conclude that radiation forces prediction can be extremely inaccurate if potential flow methods are applied to study motions of multihull configurations close to resonance frequencies. Even though recent analyses on three-dimensional problems have been presented using a potential flow based numerical method (see, for instance, Ref. [24]), in this study, we propose a numerical investigation on radiation forces of two-dimensional floating twin-hulls sections forced to oscillate at the free surface with the ultimate goal of characterizing the dependence of hydrodynamic coefficients from oscillation amplitude. A simple catamaran-like section studied by Ohkusu [25], Wang and Wahab [26], and Lee et al. [27] has been used to validate and verify the numerical model here proposed. The exact prediction of the hydrodynamic radiation forces requires a correct modeling of the energy dissipated by means of outgoing waves and dissipated through vortex shedding due to viscous separation, for this reason we used a viscous flow solver with fully nonlinear free surface capabilities for the prediction of the flow around oscillating twin sections at resonance. This method is exempt from any empiricism; hence, it is able to follow the real physics also at resonant conditions. A combination between grid resolution in the free surface region and a proper tuning of the VOF parameters allows for a correct prediction of wave breaking and overtopping that may occur in the resonance frequency range. The Navier–Stokes solver has been included in a general ship motion prediction framework initially developed and validated for single sections of generic shapes (see Ref. [28]), including small waterplane area twin hull (see Ref. [29]) and then extended to twin sections (see Ref. [30]). A significant improvement with respect to traditional potential flow theories is achieved using the proposed method for radiation forces predictions in resonant gap problems; for this reason, it is used as a virtual wave flume in this study with the ultimate goal to perform a parametric analysis particularly addressed to piston-mode motions description. In the present paper, this resonance mode is in fact deeply investigated for three different geometries and four different oscillating amplitudes. Interesting comparisons are presented aimed to characterize free surface nonlinearity due to high oscillation amplitudes and resonant free surface motions.

## Numerical Model

*p*) is here obtained using an algebraic multigrid method (GAMG (geometric algebraic multigrid)) with a diagonal incomplete Cholesky smoother. The velocity term has been solved using a PBiCG method (preconditioned biconjugate gradient [31]) with a diagonal incomplete lower-upper preconditioner. One of the main characteristics of the forced oscillation problem is the mesh morphing algorithm used to diffuse boundary motion through internal grid points. An Earth fixed reference frame is used with boundary fitted grid: a sinusoidal displacement is prescribed at the moving boundaries (hull surface) and a pseudo-solid Laplace smoothing equation distributes the motion of the moving boundaries to the internal grid points through a diffusivity coefficient

*γ*following Eqs. (1) and (2) as in Refs. [32,33]:

Equation (2) describes the deformation of an initially valid mesh ($tn\u22121$) into a deformed mesh at *t _{n}*. The quality of the grid has been checked at each time step to ensure that no faces and no cells were inverted during motion. The initially valid grid remained a valid grid, demonstrating the accuracy of the motion distribution law. The numerical solution is achieved using a preconditioned conjugate gradient with a diagonal incomplete Cholesky preconditioner. The pressure–velocity coupling is here solved using the pressure implicit with split operators scheme of Issa [34].

*α*for which a transport equation is solved together with Navier–Stokes equations

*α*close to the free surface through a relative velocity field between water and air defined as

where *n _{f}* is the face normal flux,

*S*the face area vector and $\phi $ the volume metric flux. The sharpness of the free surface is controlled through

_{f}*C*, an empirical coefficient usually varying from 0 (smeared interface) to 1 (sharp interface). In this study, a

_{f}*C*equals to 0.5 has been used to ensure a sharp outgoing wave profile prediction (obtained with high values of

_{f}*C*) and numerical stability (obtained with small

_{f}*C*values). The unsteady solution of the Navier–Stokes equations is achieved differentiating in time with an implicit Euler scheme with adjustable time step selected according to the maximum Courant number (Co < 0.5 at the free surface and Co < 1 everywhere else). Navier–Stokes equations have been solved including gravity effects but dropping the hydrostatic component from the pressure field.

_{f}*ω*integrating the solution of Navier–Stokes equations in terms of normal (pressure) and tangential (shear) stresses over the hull surface. In a linear theory framework radiation forces time signals ($fHD3$) can be decomposed in a component in phase with acceleration and a component in phase with velocity. For small amplitude motions, following the linear theory and ignoring superimposition effects the in-phase component is proportional with acceleration ($\eta 3..$) through added mass coefficient

*a*

_{33}whereas the out-phase component is proportional with velocity ($\eta 3.$) through damping coefficient

*b*

_{33}

*ξ*

_{3}and an oscillation frequency

*ω*

_{3}, added mass (

*a*

_{33}) and damping coefficients (

*b*

_{33}) are obtained as follows:

In Eq. (6), the first index indicates force direction (vertical) and the second indicates motion direction (heave). Variable $FHD3(\omega 3)$ represent the amplitude of the force component corresponding to the oscillation frequency *ω*_{3}. The phase angle between vertical force and heave motion has been indicated with $\epsilon 3$.

## Verification and Validation of the Numerical Model

The viscous flow model is validated against the experiments on twin circular sections studied by Ohkusu [35], Lee et al. [27], and Wang and Wahab [26]. In the present study, *a* is the draft of the catamaran and it corresponds to the demi-breath of each hull. Different hull spacing 2*b* have been considered: $b=1.5a,\u2009b=2a$, and $b=4a$. Main dimensions of model setup are shown in Fig. 1. Considering that the radiated wavelength (for infinite depths) is proportional to the oscillation period squared according to the dispersion equation (7), each oscillation frequency has been simulated using a different computational grid in order to keep constant grid resolution in the free surface region

This ensures a consistent accuracy level in the prediction of outgoing free surface waves as well as in the prediction of boundary layer flow and vortex shedding in the near body region. In the framework of the present study, a multiblock grid has been used. Structured quads in the near body region as well as in the free surface region have been used to increase the accuracy of the numerical solver while tetrahedral elements have been used in the far field with the goal to reduce the number of cells. Figure 2 shows a snapshot of the numerical grid.

Geometric progression is used to increase the resolution close to the body (where velocity gradients are higher) and to decrease resolution close to domain boundaries to damp out the radiated waves and limit their reflection.

A slip wall boundary condition has been implemented at the bottom wall of the computational domain positioned at a distance from the free surface equal to eight radiated wavelengths. Outlet boundaries have been obtained setting a zero gradient condition for the velocity field and a zero dynamic pressure at the top and side walls of the computational domain which are, respectively, positioned at a distance equal to four wavelengths from the free surface and 12 wavelengths from the hulls centerline. Velocity field is initialized at zero in the whole computational domain except for the hull surface, where a sinusoidal motion has been imposed. Volume of fluid scalar function has been set consistently with the free surface vertical position (i.e., *α* = 1 in the water region and *α* = 0 elsewhere). A mesh sensitivity study has been performed to select cells distribution able to ensure good accuracy reducing the computational effort. Table 1 presents results obtained at oscillation frequency $6\u2009rad/s$ and oscillation amplitude $\xi =1.27\u2009cm$. This condition corresponds to a nondimensional oscillation frequency of $(\omega 2a/g)=0.56$, which represents the piston resonance mode for this particular geometry, as demonstrated in Fig. 3. Figure 2 shows the resolution in the whole computational domain as well as an enlarged view close to the hull's surface. Results obtained in terms of nondimensional added mass and damping coefficients clearly show grid independence for the three tested configurations. The resolution of mesh M2 is able to capture free surface profile even for longer radiated waves, where M1 showed numerical instabilities mostly due to high aspect ratio cells at the boundary of the computational domains. Consequently, mesh M2 has been selected. The twin-hull section has been forced to oscillate in heave motion in a wide range of frequencies. Two different oscillation amplitudes have been imposed in the experimental procedure: $\xi =0.635\u2009cm$ has been used in the low frequency range, while at high oscillation frequencies, the motion amplitude has been increased up to $1.27\u2009cm$. The midfrequency range tests were made at both amplitudes: the linear assumption between motion and forces has been verified by comparing results obtained with two different oscillation amplitudes far from the resonance peak. In resonance condition, hydrodynamic forces nonlinearly depend on the motion amplitude. The Navier–Stokes solver validation study is presented in Fig. 3 where added mass and damping coefficients have been calculated and compared with experiments in a wide frequency range, in order to cover two resonance modes in case of the largest hull separation. Two oscillation amplitudes have been simulated. Figures 3(a)–3(f) present results for different hull separations: $b=1.5a,\u2009b=2a$ and $b=4a$. A potential flow solution based on Frank method [36], corrected for viscous effects, is plot against CFD and experimental results.

Experimental measurements obtained for smaller hull separation ($1.5a$, 2*a*) may be affected by some uncertainty, since additional aluminum angle beams were added to increase the rigidity of the test rig. Added mass and damping coefficients have been derived using a classical linear approach as from Eq. (6).

The proposed CFD computations solve the indetermination of the potential flow solver, especially close to the resonance where traditional potential flow theory usually fails to give accurate numerical solutions. The drastic under-prediction of damping resulting from potential flow methods leads in fact to the overestimation of the piston mode effect unless opportunely corrected. Lee et al. [27] used model test results to smooth potential flow forces predictions, arbitrarily imposing a threshold to the magnitude of radiation force at the first resonance frequency (corresponding to the *piston mode*). Higher-order modes (*sloshing modes*) were not corrected. The proposed Navier–Stokes solver is proven to give particularly robust predictions at the resonance peak, correctly showing the effect of motion amplitude and in many cases getting close to the experimental measurements. In the low frequency range, both CFD and potential flow methods present discrepancies with experiments: Wang and Wahab [26] attributed the deviation of the experimental results to the existence of a disturbance caused by end plates lateral oscillations.

## Piston Mode and Higher-Modes Resonance Frequencies

A resonant vertically oscillating water column is experienced in vertical gaps defined by the inner walls of twin hull configurations. Peaks in the damping coefficient and abrupt change of phase of the internal wave amplitude are experienced at particular frequencies which have been closely investigated in the present study. Molin [7] defined two types of resonant motions of the inner fluid:

*Sloshing mode*: a standing wave of length proportional to hull spacing is generated between the inner sides of the catamaran.*Piston mode*: the free surface in between the inner sides of the catamaran oscillates in the vertical direction, resembling a solid body of water (piston).

In the first case water motion inside the gap presents analogy with fluid sloshing in closed tank: water behaves like there is no communication between inner and outer side of the hull section, while in the second case, mass fluxes between the inner water column and the water outside the gap are responsible of a pistonlike motion of the free surface.

In Eq. (8), *b* is the inner side distance (i.e., gap width), while *a* is the draft, $B=2(b+a)$ is the overall breadth of the barge, and $H=\lambda B$ (with $\lambda >1$) represents the position of the sinks used to simulate the water flow into the moonpool.

In the numerical method presented in this paper, the free surface solution is achieved through a surface capturing technique, so nonlinear effects, for instance, wave breaking, are inherently included in the flow solution, contrary to the method presented by Molin [7] that uses linearized boundary conditions on the free surface. Figure 4 demonstrates the importance of vortex shedding and viscous effects: in each snapshot, the $b=1.5a$ catamaran section has been forced to oscillate at four different amplitudes. Generally, the amplitude of the free surface oscillation inside the two hulls increases with the motion (from top to bottom), until it breaks. In Fig. 4, the vorticity field has been nondimensionalized multiplying dimensional values by $a/\xi \omega $. Fully potential-flow methods or potential-flow methods combined with vortex-tracking strategies at the bottom opening cannot properly consider the interaction between the vorticity shed by hull motion and the free surface if a viscous flow solver is not included in the physical model. Nonlinear effects due to piston mode resonance of the free surface inside the two hulls are evidenced by the coalescence of vorticity at the free surface as shown in Fig. 4. Alternate vortex shedding is persistent in the fluid domain and captured in Fig. 4(d) where multiple vortices detaching from the keel and convected toward the outer free surface are visible. No symmetry conditions have been implemented at the catamaran symmetry plane.

Nonlinear behavior of piston resonance mode has been studied considering a systematic variation of the forced heave amplitude (*ξ*) and hull spacing (*b*) over a wide frequency range, characterizing the dependence of radiation forces (added mass and damping coefficients) on the investigated parameters. At each time step, the free surface has been identified through the vertical position of the cell centroid having $\alpha =0.5$ : four different probes (A, B, C, and D in Fig. 1) have been used to measure the wave elevation of the internal (probe A) and the outgoing waves at different distances from the hull (B, C, and D).

Figure 5 shows an example of the calculated free surface elevation versus time at the four probes locations at an oscillation frequency $\omega a/g=0.81$. Free surface time histories presented in Fig. 5 correspond to an oscillation frequency close to resonance as evident from the comparison of the internal wave amplitude (A) and the radiated waves (B, C, and D): internal wave amplitude *ζ* is in fact about 10 times larger than the motion amplitude *ξ*.

Figure 6 presents at a glance the most significant results of the parametric investigation: the whole figure is composed by 12 graphs arranged in a 4 × 3 matrix. Each row corresponds to a given heave amplitude, while each column corresponds to a fixed hull distance (gap width). In each plot, the nondimensional added mass coefficient (continuous line) and the nondimensional damping coefficient (dashed line) are presented together with the response amplitude operation (RAO), i.e., the ratio between the internal wave amplitude *ζ* and the motion amplitude *ξ* as a function of frequency (dashed–dotted line). The nondimensional frequency range goes from $\omega a/g=0.249$ to $\omega a/g=1.496$. A sufficiently small frequency step ($\Delta \omega =0.25\u2009rad/s$) has been selected in order to obtain smooth variations, especially where strong gradients due to resonance effects were expected.

The general hydrodynamic characteristics and nonlinear nature of piston or sloshing resonant modes can be inferred from Fig. 6. The rapid change of phase between motions and radiation forces which causes the added mass coefficient to assume negative values corresponds to an inversion of the trend of the damping coefficient which shows, at this frequency, a maximum value. For slightly higher frequencies, a very large inner wave amplitude is experienced (first peak of the dashed–dotted line). If the frequency is further increased, the damping coefficient rapidly decreases featuring a minimum value corresponding to the so-called *trapped mode*. Figure 6 demonstrates the existence of resonant phenomena even for larger hull distances.

Sloshing resonance modes are experienced also at lower hull separations at frequencies out of the investigated range. Piston and sloshing modes are a consequence of two distinct phenomena: the former corresponds to a maximum mass flux between the inner and the outer fluid regions while the latter is experienced when the flow exchange reaches its minimum value. Piston mode resonant condition is experienced around a single low frequency, while sloshing modes are found at discrete number of higher frequencies proportional to hull distance.

It is interesting to compare the first resonant frequency found from the nonlinear viscous model with the theoretical values obtained by Molin for a twin rectangular section, which are often used for practical purposes. These frequencies are obtained applying Eq. (8) originally derived for moonpools to catamaran geometries. Values in Table 2 are reported in Fig. 6 by means of dashed vertical lines for the three different hull spacings. Differences in hull shape have to be considered when comparing values obtained for a moonpool section (Eq. (8)) and the one predicted for twin circular cylinders.

Figures 7(a) and 7(b) present results at a fixed oscillation amplitude $\xi =1.27\u2009cm$ : nondimensional added mass and damping coefficients are plotted for the three different hull spacings together with results obtained for a single hull. Numerical predictions presented in Fig. 7 show consistent limit trends with the expected physics: in the high frequency limit, the nondimensional added mass coefficients converge to 1 while the nondimensional damping coefficients converge to 0 (same of the single hull case).

Figure 8 presents results for $b=1.5a$ oscillating at $\xi =1.270\u2009cm$ in terms of added mass and damping as well as the nondimensional internal wave elevation at the four probes previously described (A, B, C, and D). This plot evidences that the maximum value for the damping coefficient occurs at the zero crossing of the added mass coefficient: radiation forces are perfectly in phase with velocity; hence, energy dissipation reaches its maximum value; this condition anticipates a rapid drop of the added mass coefficient. Slightly increasing the oscillation frequency leads to the trapped mode: a persistent fluid motion is experienced in between the hulls. Energy is not dissipated through wave generation since no waves are radiated as demonstrated by the vanishing outgoing wave amplitude probed at B, C, and D presented in Fig. 8. For $b=1.5a$, the zero damping condition is experienced at $\omega =7.25\u2009rad/s$ and corresponds to a nondimensional frequency of $\omega a/g=0.903$. For this oscillation frequency, the internal flow motion interacts with the radiated flow leading to an almost perfect outgoing wave cancelation. In this condition the fluid volume displaced by the hulls during their motion nearly equates the volume change of the water column in between the two hulls. Since energy losses due to wave generation are the most important source of damping and being energy dissipation significantly dependent on outgoing wave amplitudes, damping coefficient falls down toward zero when this particular frequency is approached and viscous effects importance becomes more relevant. In Fig. 9, volume of fluid fraction contours and wave elevation plots are shown for two different oscillation frequencies, corresponding to the piston mode frequency ($\omega =6.25\u2009rad/s,\u2009\omega a/g=0.779$) and to the zero damping frequency ($\omega =7.25\u2009rad/s,\omega a/g=0.903$).

The piston-mode motion phenomenon is analogous to the problem of an oscillating water column in between the two hulls. We present the solution to this problem with the specific aim to confirm and verify the ability of the proposed numerical model to capture the characteristics of resonant flows dominated by nonlinear viscous effects. The height of the water column has been initially set to the maximum level numerically predicted at the piston mode resonance condition and let free to oscillate in order to find a stable condition as in a decay test. Wave elevation has been recorded at each time step using a probe located at position A. The time history of the water level in Fig. 10(a) clearly evidences the oscillatory behavior of the water column. The amplitude decay is mainly due to wave damping. The analysis of the periodic irregular time history of the water column between the two hulls has been performed using a FFT algorithm as presented in Fig. 10(b). The principal component of the FFT is found at $\omega a/g=0.803$. This value has been compared with those found with forced oscillation of the same catamaran configuration ($b=1.5a$ in Fig. 11). A very good match between the characteristic frequencies of this phenomena has been found at smaller oscillation amplitudes: for $\xi =0.635\u2009cm$ the resonance frequency is $\omega pisa/g=0.807$, while for $\xi =1.27\u2009cm$ it is $\omega pisa/g=0.792$. The correspondence with the higher amplitudes is weaker because of free surface rupture (as seen in Figs. 4(c) and 4(d)).

### Oscillation Amplitude Effects.

The extensive simulation campaign presented in Fig. 6 allows to investigate over of the inner wave amplitude and the piston mode frequency dependency on motion amplitude for different hull separations. Results presented in Figs. 11 and 12 show for different oscillation amplitudes *ξ* the trend of piston mode frequency and nondimensional inner wave amplitude (RAO), respectively. A consequence of nonlinearity is demonstrated in Fig. 11: for shorter hull distance piston mode frequency is significantly dependent on motion amplitude, in particular for larger oscillation amplitudes resonance is experienced at shorter oscillation frequencies. This nonlinearity decreases for larger hull separations until it almost disappears for $b=4a$ for which an almost constant resonance frequency has been found. Another evident result is that piston mode is experienced at shorter frequencies for larger hull separations. Figure 12 shows results in terms of maximum water column height as function of heave motion amplitude (RAO): one curve has been presented for each hull configuration. The nondimensional maximum inner water level has a convergent trend for high oscillation amplitude. Moreover, when hull interference increases (smaller gap widths), the RAO is enhanced at a lower oscillation amplitude; this is evident for $b=1.5a$ and $b=2a$ configurations (continuous and dashed line in Fig. 12). An opposite trend has been found for larger hull separations ($b=4a$) for which RAO is slightly higher at higher oscillation amplitudes.

Resonance effects due to hull interference are more significant for closer hull configuration; in these cases (e.g., $b=1.5a$), the RAO is significantly reduced at higher oscillation amplitudes. At lower oscillation amplitudes, energy dissipation due to viscous effects is in fact lower; hence, fluid motion in between the hulls is not significantly damped out. Moreover, higher oscillation amplitudes lead to higher flow velocities responsible for a nonlinear dynamic of the free surface which eventually results in wave breaking (e.g., Figs. 4(c) and 4(d)). At lower amplitudes, the internal water column heaves with a lower velocity, hence reducing wave breaking occurrence. Resonance effects and nonlinearity become weaker for larger hull separation ($b=4a$) for which wave breaking has not been experienced at any tested oscillation amplitude: for this configuration, the RAO follows a different trend, actually increasing with the motion amplitude (dashed–dotted curve in Fig. 12). This demonstrates weaker piston mode effects, actually characterized by higher horizontal flow components.

Figure 13 presents the RAO at piston mode frequency with respect to hull distance *b*/*a*: one curve for each oscillation amplitude is presented. Same conclusions can be drawn looking at Fig. 13: when the gap width becomes small (*b* = 1.5*a*), the RAO significantly depends on the forced oscillation amplitude while for larger gaps, the RAO asymptotically converges to unity. In the piston mode frequency range, the fluid behaves like a solid body; hence, the forced heave oscillation causes part of the displaced fluid to be moved outside and part inside the hulls. Figure 13 evidences that for smaller gap widths, a higher water column is expected in between the two hulls. If the overall trend of all curves is taken into account, it is evident that increasing gap width (*b*) leads the inner wave amplitude to convergence to the same value regardless of motion amplitude.

## Conclusions

In this study, we address resonant nonlinear flow motion around twin oscillating sections by a fully viscous free surface flow solver.

A viscous numerical solver based on the solution of the unsteady Navier–Stokes equations has been presented together with a surface capturing technique based on the solution of the VOF scalar function. Both outgoing and internal waves must be predicted with very high fidelity in order to obtain good accuracy even for very large free surface deformations which eventually lead to wave breaking. The VOF approach allows for the solution of Navier–Stokes equations for a single fluid mixture: no boundary conditions are required at the free surface; hence, nonlinearities can be directly included in wave dynamic predictions. The numerical method has been first validated against Wang and Wahab [26] experimental results on a catamaranlike section often used as a benchmark study for ship motions prediction codes.

Particular attention has been focused in resonance frequencies where even for small amplitudes the linearized free surface condition often used in potential flow methods (current state of the art seakeeping codes) represents a strong approximation. The numerical method is able to solve equally well the two different resonant modes: the piston mode motion occurring at a single lower frequency corresponds to the first resonant mode while the sloshing modes motion correspond to higher modes and occur at frequency multiple of hull distance.

Traditional potential flow theories are not suitable to solve motion amplitude and radiation forces around resonance, whereas the viscous method with nonlinear free surface lead to very accurate results even for first mode resonant conditions where the jet-flow type induced by piston-mode effect causes very high water flow motions.

At this particular frequency, there is in fact a significant phase change of hydrodynamic forces: the added mass rapidly becomes negative while the damping coefficient reaches its maximum value. Just beyond this frequency, the outgoing wave vanishes and so does the damping coefficient.

A very significant dependence of added mass and damping coefficient on the amplitude of the forced oscillation has been found close to resonance; for this reason, an extensive simulations campaign has been performed with the ultimate goal of characterizing hydrodynamic coefficients dependence on oscillation amplitudes. Almost 500 two-dimensional simulations have been performed: four different heave amplitudes have been considered in order to characterize unsteady hydrodynamic forces around oscillating twin sections with three different spacings. Nonlinear correlations have been found. A first straightforward conclusion is that radiation forces dependence on oscillation amplitudes is particularly enhanced by hull proximity, including viscous effects and nonlinear treatment of the free surface allows on one hand to predict vortex shedding and its interaction with the free surface and on the other hand to capture wave breaking occurring at higher oscillation amplitudes. The presented parametric analysis demonstrates that high relative piston mode motions are found at lower motion amplitudes. Moreover, it has been found that first mode resonance frequency occurs at shorter oscillation frequency for larger oscillation amplitudes.

Finally, the method has been validated also through a comparison with the natural decay frequency of an oscillating water column initialized with an unbalanced level between inner and outer free surface. Also in this case, from numerical simulations we were able to derive the exact natural frequency of the oscillating fluid motion between the hull, which is known to be very dependent on viscous damping. The consistency of the numerical method presented in this paper is confirmed by the exact matching of the natural frequency for these two analogous problems.

This paper ultimately demonstrates that piston mode frequency and internal wave height strongly depend on the motion amplitude and on the gap width in between the two hulls. As a final conclusion, the viscous flow method with fully nonlinear free surface leads to higher fidelity results in case of resonant fluid motions in gap problems and it is able to capture the nonlinearity induced by high amplitude motions.

## Acknowledgment

The authors wish to express their gratitude to Dr. Fariba Fahroo for the support received by the DARPA EQUiPS grant.

## Funding Data

Defense Advanced Research Projects Agency (Grant No. HR0011517798).