Coupled Effects of Incident Waves Forcing and Internal Tank Sloshing on the Dynamics of Twin Floating Bodies

1 Introduction

The problems of tank sloshing and multi-floating bodies hydrodynamics are research hotspots in the field of naval architecture and ocean engineering for both academic research and engineering applications. For example, offloading of LNG can be undertaken in a tandem arrangement, where an LNG carrier is moored at the stern of the FLNG vessel [1], or in a side-by-side configuration where the LNG carrier is moored alongside the FLNG vessel. When the moored floating bodies are in restricted water areas, resonance responses of the gap flow or multi-bodies may be induced at specific wave conditions, which can result in a significant increase in the wave dynamic effect or wave height in narrow gap areas. The increase in wave dynamic effect can also increase the motions and wave loads acting on the floating bodies in return. This will pose a threat to the structure strength, collision and operation safety. Therefore, the hydrodynamic disturbance of the multi-floating body structures has received increasing attentions [2].

The research on the hydrodynamics of water wave interaction with multi-bodies has been widely conducted. For example, Saitoh et al. [3] experimentally investigated the fluid resonance effect in the gap between two fixed bodies in a 2D plane. Fredriksen et al. [4] analyzed the dynamic variation of wave height in a narrow gap between two moving boxes by experimental method. Zhu et al. [5] investigated the interaction force of two floating bodies with different gap widths by potential flow theory. However, the potential flow theory tends to overestimate the resonance amplitude at the gap resonant frequency, which is due to the neglect of fluid viscosity and the nonlinear free surface effect. To address this problem, Tan et al. [6] proposed a theoretical dynamic model that includes a linearized damping coefficient, and the obtained numerical result agrees well with that by viscous fluid simulation. Zhang et al. [7] studied the effect of gap width on the flow behavior of 3D fixed double floating bodies using finite difference method of viscous flow theory. Lu et al. [8] investigated the resonance responses in an array of three bodies, and the responses at the two gaps are studied, which is significantly different from the performance in a single gap of twin floating bodies configuration. Similar work can also be found in Sun et al. [9], Jing et al. [10] and Gao et al. [11]. However, in these studies, the multi-bodies are fixed structures, where the radiation problem is not considered, as the wave-induced movement of the structure is not involved.

In fact, many ships and offshore structures are floating structures, and the wave-induced movement of structures affects the flow field in return. Therefore, the interaction of multi-floating bodies and the gap flow resonance is also a reason for concern. He et al. [12] studied the wave-induced fluid resonance in a narrow gap between two rigid-connected boxes that have only one degree of freedom of heave released. Jing et al. [13] further investigated the effects of heave and roll motions on the resonant frequency and amplitude of the gap flow. Sun et al. [14] studied wave-induced free surface motion in the gap between a tanker and an FLNG barge under both fixed and free-floating conditions. Hong et al. [15] experimentally studied the hydrodynamic behavior of Tandem and side-by-side moored vessels of LNG carrier and FLNG in different wave headings. Fang and Chen [16] studied the relative motion and wave elevation in the gap between two ships advancing in waves, and the calculated results are well validated by experimental data. Choi and Hong [17] used a high-order boundary element method to accurately calculate the hydrodynamic interaction between multiple floating bodies, which includes ship motion response and wave drift force. Hong et al. [18] studied the hydrodynamic interaction of side-by-side moored multiple vessels by both numerical and experimental methods. More studies about the gap resonance of double floating bodies can be found in Jing et al. [19], He et al. [20] and Jing et al. [21]. Mohapatra et al. [22] presented a numerical hydrodynamic model for a moored interconnected floating large structure under the action of regular waves to analyze the effect of current and wind. A review of recent developments on the hydroelastic response and gap resonance of multi-body floating structures is presented in Amouzadrad et al. [23].

Moreover, sloshing in the liquid tanks makes the liquid cargo ships’ motion responses much different from the solid loading condition. To date, the coupled response of tank sloshing and floating body motions has been studied by different methods, including potential flow theory [24], CFD simulations [25], SPH simulations [26] and model experiments [27]. Although multi-floating body systems and liquid tank sloshing have been widely investigated, these two problems are usually studied separately. The fully coupled analysis of liquid tank sloshing inside floating bodies in side-by-side configuration is rarely conducted. Jiao et al. [28] simulated two side-by-side LNG ships’ motions coupled with tank sloshing in regular waves by SPH solver. Zhao et al. [29] developed a potential flow theory code to predict the dynamic responses of side-by-side arranged FLNG system, and solved the coupled responses between sloshing in liquid tanks and vessels’ motions in the time domain.

Currently, research on the hydrodynamic interaction problem of twin-floaters with liquid tanks remains quite limited, especially regarding the influence of liquid sloshing in the tanks on the gap flow within the twin-floaters system under multi-degree-of-freedom motions, which awaits further exploration. To better understand the resonance mechanism and response characteristics of twin-floaters under regular waves, this paper conducts a comprehensive study on the hydrodynamic behaviour of simplified floating body of both one- and twin-floater configuration and also with and without tank sloshing. The paper is organized as follows. The basic method of SPH algorithm is described in Section 2, and the applicability of the SPH method in solving the twin-body hydrodynamic problem is also validated by comparing with reference data. In Section 3 the numerical model setup of the twin-floater with tank sloshing involved in this paper is established, and the convergence analysis of time step and particle spacing is conducted. Comparison of responses between one- and twin-floater systems for the cases of without and with internal liquid tank is presented in Section 4. In Section 5, the results are comprehensively compared and analyzed.

2 Basic Theory and Algorithm Validation

This study uses the DualSPHysics open-source software platform, which is on the basis of smooth particle hydrodynamics (SPH) to implement the hydrodynamic simulations. The simulation calculation of the multiphase fluid flow problem is performed via the SPH algorithm of graphics processing unit (GPU) acceleration technology. The basic theory of SPH algorithm is briefly described first. The wave generation and absorption method in the numerical wave tank using the SPH tool is presented. Then the hydrodynamic problems of two fixed boxes and the gap flow resonance phenomenon are simulated by the developed numerical wave tank, and the results are validated by comparing with experimental data in reference.

2.1 SPH Theory

Using the integral equation on the basis of interpolation function, the conservation law of continuum fluid dynamics is converted into the partial differential form suitable for particle simulation, and an estimated value at a specific point is given. This interpolation or weighting function is usually called the kernel function W, which can be expressed in different forms. The most commonly used form is the cubic or quintic function. It is generally expressed as a function of F(r) defined by r′: F(r)=∫F(r′)W(r−r′,l)dr′(1) where r is the displacement vector of the space point, r′ is the displacement vector of the adjacent space point, and l is the smoothing length, W(r − r′, l) is a smoothing kernel function which is an equivalent function that approximates the properties of the Dirac function. Fig. 1 shows the center of particle a and the radius 2l of the integration region involved in the calculation.

Figure 1: Sketch of the particle integration region.

The performance of the SPH model largely depends on the selection of the smoothing kernel function. In this work, the cubic spline kernel function used in Wendland [30] is as follows: W(r−r′,l)=αD(1−q2)4(2q+1) (0≤q≤2)(2) where q = r − r ′ l and α D is a coefficient in the 2D calculation α D = 7 4 π l 2 ; and in the 3D calculation α D = 21 16 π l 3 .

The equation of momentum conservation in the continuum is as follows: dvdt=−1ρ∇P+g+Γ(3) where v represents the fluid velocity, ρ represents the fluid density, P represents the normal pressure, g represents the acceleration of gravity, and Γ represents the dissipation term. DualSPHysics provides more options for expressing the dissipation term of artificial viscosity. The scheme has certain simplicity, and it is a common method for fluid simulation via SPH. In this work, the artificial viscosity scheme (Monaghan [31]) is used to express the dissipation term, and the momentum conservation equation can be expressed by the following equation: dvadt=−∑bmb(Pb+Paρb⋅ρa+Πab)∇aWab+g(4) where P a , P b are the pressure of particle a and b, respectively; ρ a , ρ b are the density of particle a and b, respectively; the viscosity term ∏ a b is expressed by the following formula: Πab=−ac¯abμabρ¯abvab⋅rab<00vab⋅rab>0(5) where v a b = v a − v b , r a b = r a − r b , and v a , v b and r a , r b are the speed and position of particles a and b, respectively. The smaller the particle spacing, the more particles are required, and the longer the calculation time will be. The computational complexity of SPH is theoretically approximately linear, but in practice, it is significantly affected by neighbor search, physical models, and hardware limitations. μ a b and c ¯ a b are represented by the following equations [32]: μab=lvab⋅rabrab2+η2(6) cab¯=0.5(ca+cb)(7) η2=αl2(8) where c a , c b are the sound speeds of particles a and b. α is an adjustable coefficient used to introduce appropriate dissipation, and α = 0.01 is used. During the entire simulation process of the SPH method, the mass of each particle remains constant; only their relative density fluctuates. The changes in the relative density of these particles are calculated by mass conservation or the continuity equation which can be expressed via the SPH method: dρadt=ρa∑bmbρbvab⋅∇aWab(9) where m b is the mass of adjacent particle b. In this work, according to the studies of Monaghan [33] and Batchelor [34], the relationship between pressure and density is expressed as: P=b(ρρ0)γ−1(10) b=c02γρ0(11) c0=c(ρ0)=(∂P/∂ρ)ρ0(12) where P is the pressure, ρ is the density, in the fluid calculation γ = 7, ρ 0 = 1000 kg/m3 is the reference density of the fluid, and c 0 is the speed of sound at the reference density. In this work, the artificial viscosity model described in Monaghan [31] is used to simulate the motion of a floating body in a numerical tank to address the viscous dissipation of momentum. The Symplectic scheme in Ren et al. [35] is used in the time-stepping algorithm, which is of second-order accurate in time and more conducive to wave generation.

2.2 Wave Generation

2.2.1 First Order Wave Generation

The Biesel transfer functions express the relation between wave amplitude and wavemaker displacement, under the assumption of irrotational and incompressible fluid and constant pressure at the free surface in Altomare et al. [36]. The transfer function links the displacement of the piston-type wavemaker to the water surface elevation, under the hypothesis of monochromatic sinusoidal waves in one dimension in the x-direction: ηx,t=h2coswt−kx+δ(13) where h is the wave height, d is the water depth, x is distance and δ is the initial phase. The quantity ω = 2π/T is the angular frequency and k = 2π/ λ is the wave number with T equal to the wave period and λ the wave length. The initial phase δ is given by a random number between 0 and 2π. This equation expresses the surface elevation at infinity. The Biesel function can be derived for the far-field solution and for a piston type wavemaker as Jiao et al. [37]: hS0=2sinh2kdsinhkdcoshkd+kd(14) where S0 is the piston stroke. Once the piston stroke is defined, the time series of the piston movement is given by:

2.2.2 Second-Order Wave Generation of Regular Waves

The Symplectic scheme in Ren and Khayyer [35] is used in the time-stepping algorithm, which is second-order accurate in time and more conducive to wave generation. The implementation of a second-order wavemaker theory will prevent the generation of spurious secondary waves. The second-order wave generation theory implemented in DualSPHysics is based on Madsen [38], who developed a simple second-order wavemaker theory to generate long second-order Stokes waves that would not change shape as they propagated. The theory proposed by Madsen is straightforward, controllable, computationally inexpensive with efficient results, and is accurate for waves of first and second order.

The piston stroke S0 can be redefined from Eq. (14) as S0 = h/m1 where:

Following Madsen [38], to generate a wave of second order, an extra term must be added to Eq. (15). This term, for piston-type wavemaker, is equal to:

Therefore, the piston displacement, for regular waves, is the summation of Eqs. (15) and (17):

2.2.3 Passive Wave Absorption

A damping zone is implemented in DualSPHysics as passive absorption system. The implemented damping system consists in gradually reducing the velocity of the particles at each time step according to their location, but using quadratic decay rather than exponential. In this way, the velocity is modified following: v=v0⋅f(x,Δt)(19) where v0 is the initial velocity of the particle i, v is the final velocity and f(x, Δt) is the reduction function defined as: f(x,Δt)=1−Δt⋅β⋅x−x0x1−x02(20) where Δt is the duration of the last time step, x is position of the particles, x0 and x1 are the initial and final position of the damping zone, respectively. It is recommended to use one wavelength, L, as the length of the damping zone. The coefficient β modifies the reduction function that is applied to the velocity.

2.3 Numerical Method Validation

To validate the SPH numerical method in the simulation of twin-floater hydrodynamics, a numerical model was established that is the same as the experimental model in Saitoh and Miao [3] (see Fig. 2). The 2D numerical tank is 16 m long and 1 m high, and the water depth H = 0.5 m. Regular waves with a predefined wave height h and wave number k are generated by a piston-type maker which is located at the left side of the tank. A damping zone is set behind the floating body at the other end of the tank. For the boundary conditions, an improved adaptive dynamic boundary condition (mDBC) is used on the wall and surface of the numerical tank to avoid defects such as excessive dissipation that may be caused by the traditional dynamic boundary condition (DBC) [39]. The SPH method uses particles to discretely distribute the fluid in the water tank. The push plate moves according to the originally set program. The water particles move accordingly under the action of the push plate, which can visually depict the wave shape at each moment. The entire simulation process is relatively in line with actual cognition and avoids the problems of grid deformation and reconstruction in the wave simulation process of the traditional grid method.

Two fixed boxing bodies are adopted so as to be in consistent with the experimental model in reference. The two fixed bodies have the same geometric size of length L = 0.5 m, depth D = 0.5 m, and draft d = 0.252 m. The width of the narrow gap between the two fixed bodies is set at s = 0.05 m. The fixed body is placed at a distance of 5 m from the wave-maker. Incident waves at different wave number k = 1.2 to 2.4 with a constant wave height h = 0.024 m are generated for study. Two wave probes are used to monitor the waves: Probe 1 is installed at X = 3 m from the wave-maker to monitor the incident wave without disturbance, Probe 2 is placed at X = 5.53 m, which is used to measure the gap flow between the two bodies. The particle spacing of 0.001 m is used in this case.

Figure 2: Schematic diagram of the construction of the numerical tank used in the numerical validation.

The case of wave number k = 1.5 and wave height h = 0.2 m is used to illustrate the generated incident waves. Comparison of the time series of wave elevation at Probe 1 in the numerical tank between the monitored and the theoretically expected values is shown in Fig. 3. It shows high consistency between the two curves, which confirms the effectiveness of the numerical wave-making method in this study. The wave height of SPH simulation is 0.11 m, which is 10% larger than the theoretical value of 0.1 m. The wave elevation by SPH simulation shows obvious nonlinear characteristics of sharp crest and flat trough compared with the linear theory wave.

Figure 3: Comparison of time series of wave elevation at Probe 1 (wave number k = 1.5).

The narrow gap flow between the two bodies reveals complex characteristics due to the resonance responses. Fig. 4 shows the ratio of the wave height hD in the narrow gap measured by Probe 2 to the incident wave height h at Probe 1 for the cases of different wave numbers. The results are compared between the present SPH method in this study, the theoretical calculation of Saitoh and Miao [3] and the experimental measurement of Lu and Cheng [8]. Notably, all the results exhibit a significant increase in wave height at narrow gap when the wave number is approximately 1.6, which indicates that a strong resonance phenomenon occurred at or near this wave number. The SPH model used in this study accurately captures this resonance effect, although the obtained wave height at the resonant position is slightly lower than the experimental value (with an error of −6%). The error may be caused by the factor of particle spacing or total number of particles as more accurate flow characteristics can be simulated with smaller particle spacing [37]. When the wave number k tends to 1, the ratio of hD/h tends to 1, which means that long waves have strong transmissivity, and the gap resonance effect becomes weak. However, when the wave number k exceeds 2.5, the ratio of hD/h tends to 0, which means that short waves have weak diffraction and the waves are reflected from the structure body.

Figure 4: Variations in the wave height in the gap at different wave numbers compared between different methods.

The time series of wave elevation at Probe 2 in the narrow gap under the resonance condition (wave number k = 1.5) is shown in Fig. 5, where the time series of the monitored incident wave at Probe 1 and the theoretical wave amplitude are also presented for comparison. It is seen that the wave height at the narrow gap reaches four times of the incident wave height, which indicates that severe resonant response occurs in the gap flow. The phase of the time series of wave elevation at the two positions is consistent. To further elucidate the mechanism of the resonance, the dynamic variation characteristics of the velocity field in the gap region within one wave period are shown in Fig. 6 for the case of wave number k = 1.5. The fluid velocity in the gap area is greater than the wave field at other area. The wave energy to concentrate on the gap region, thereby promoting an increase in fluid velocity. This increase in the velocity field leads to an increase in the wave height in this region.

Figure 5: Time series of the wave elevation in the gap and the undistributed incident waves under the resonance wave excitation case (wave number k = 1.5).

Figure 6: Velocity field distribution in the gap under resonant wave excitation case (T₀ = 10 s).

3 Numerical Model Setup of Twin-Floater with Tank

3.1 Model Description and Boundary Conditions

To investigate the hydrodynamic behavior of 2D floating body in regular waves, a 2D numerical wave tank was established, which is shown in Fig. 7. To study the hydrodynamic characteristics of different floating body models, four floating body schemes are adopted: single floater (Scheme 1), twin-floater (Scheme 2), single floater with liquid tank (Scheme 3) and twin-floater with liquid tank (Scheme 4). All the floating bodies in different schemes have the same dimension of length L = 0.5 m by depth B = 0.2 m with draught of T = 0.1 m. For the floating body without liquid tank, a uniform solid model is used, which has an average density of 500 kg/m3 and a pitch moment of inertia of 1.04 kg·m2. For the floating body with liquid tank, it has a cavity structure inside with a wall thickness of 10 mm. The density of the wall structure is 500 kg/m3. The volume of the loaded liquid is half of the whole cavity volume. The fluid inside the tank is water, which has a same density of 1000 kg/m3 with the external fluid.

Figure 7: Sketch of the numerical wave tank: (a) overview of the numerical tank and (b) floating body models.

The distance between the center of gravity of Floater 1 and the wave-maker is X = 5 m. For the cases of the twin-floater model, the gap distance between the two floating bodies is s = 0.2 m. For the cases of the single floater model, its location coincides with the Float 1, which is the front side one of the twin-floater. Only two degrees of freedom of heave and pitch are released for all the floating bodies. Two axial forces of the horizontal and vertical directions are monitored. Fig. 8 shows an illustration of the two degrees of freedom and force directions that are defined in the local coordinate system G-xz. The horizontal force represents the fluid force or fixing force acting on the body, while the vertical force represents the total force of fluid and gravity on the body.

Figure 8: Definition of the movement and force of the floating body. (a) motion; (b) force.

To better extend the generality of the conclusions, the force is dimensionlessized by: Fx¯=FxρghDL(21) Fz¯=FzρghDL(22) where Fx and Fz are the horizontal and vertical forces on the floating body, respectively; ρ is water density, g is gravitational acceleration, h is wave height, D is the body’s depth, L is length of float body.

3.2 Calculation Conditions

In this study, regular waves with constant wave height of h = 0.1 m are used, while the response of the floating body under different cases are involved by changing the wave length or wave number k. Table 1 lists the detailed calculation conditions. Table 2 explains the corresponding wavelengths and wave periods of each wave number.

For all the calculations, an AMD3990X 64-core CPU processor and an RTX3090 24 GB GPU were used. The number of particles is set to 488,785. The simulation time step size is 0.05 s. The calculation time consumed for each condition is approximately 1.5 h.

3.3 Verification and Convergence Analysis

Both the time step size and the particle spacing affect the convergence of the results. For the particle spacing dp, three sets of spacings are selected: dp = 0.01 m, 0.004 m and 0.002 m. Fig. 9 compares the time series of motions of Floater 1 of the twin-floater system (Scheme 3) that obtained by the three particle spacing schemes for Case 16 with wave number k = 1.7 and wave height h = 0.1 m. The motion state of the floating bodies and flow distribution by the three different particle spacings at a same time instant are compared in Fig. 10. This reveals that the result by the medium particle spacing dp = 0.004 m agree well with the fine particle spacing dp = 0.002 m. Therefore, the convergence analysis on the time step size is further conducted using the medium particle spacing dp = 0.004 m.

Figure 9: Convergence verification of particle spacing (Case 16).

Figure 10: Motion state of the floating bodies and flow distribution by different particle spacings (Case 16).

The motion response results by three sets of time step, i.e., 0.05 s, 0.01 s, and 0.005 s, are compared in Fig. 11. The results indicate that the calculated motions are highly consistent between the three different time step sizes. As a compromise of the calculation efficiency and accuracy, the parameters of particle spacing dp = 0.004 m and time step size dt = 0.05 s are chosen for the subsequent numerical calculations.

Figure 11: Time step convergence verification (Case 16).

4 Analysis and Comparison of the Results

In this section, the hydrodynamic behaviours of single- and twin-floating bodies are comparatively analyzed, which mainly concentrate on the heave and pitch motions as well as the horizontal and vertical forces on the floating bodies. The analyses are conducted for the cases of floaters without and with liquid tank.

4.1 Comparison between One- and Twin-Floater without Liquid Tank

In this part, the flow characteristics and dynamic behavior of the free-floating bodies under wave excitation are investigated. The response is compared between the one- and twin-floater systems.

4.1.1 Motion Responses

The variation in the motion amplitude of the floating bodies without liquid tank under different incident wave cases is shown in Fig. 12, where the results of the three floaters are presented for comparison. In general, for all three floaters the largest motion response amplitude is at k = 1.55 (corresponding to period T = 1.18 s) where the resonance responses of the floating body system occur. The resonance period is close to the natural period of the water in the narrow gap Tg = 1.11 s which is obtained according to the study of Jing et al. [21]. The motion response of Floater 2 in the twin-floater system is the lowest among the three floaters over the whole wave number range. This is caused by the incident wave shielding effects due to the presence of the Floater 1, which is located in front of the Floater 2 [40]. The difference in the heave motion of Floater 1 between the one- and twin-floater systems is small. However, the pitch motion of Floater 1 in the twin-floater system is much greater than that of the single one, especially at the resonance wave range. As the wave number increases, the pitch amplitudes of Float 1 are close for the one- and twin-floater systems.

Figure 12: Motion amplitudes of the floating bodies without tank at different wave numbers.

Fig. 13 and Fig. 14 compare the time series of motion between the three floaters for the two cases with the wave number k = 1.55 and k = 2.2, respectively. It is seen that there is an obvious phase lag in the motion signal between the Floater 1 and Floater 2 in the twin-floater system, while the phase is almost the same for the Floater 1 of the one- and twin-floater systems. This shows that in a twin-floater system, the gap flow has an important effect on the floating body at the rear (Floater 2), but has little effect on the floating body at the front (Floater 1). The motion amplitude value of Floater 2 is much lower than that of Floater 1 in the case of k = 2.2, which is also reflected in the amplitude variation results. This shows that both the amplitude and phase of the rear floating body (Floater 2) in the twin-floater system are affected by the gap flow.

Figure 13: Comparison of the motion time series of the floating bodies without tank in the case of k = 1.55.

Figure 14: Comparison of the motion time series of the floating bodies without tank in the case of k = 2.2.

To provide a deeper understanding of the resonance mechanism, Fig. 15 presents the distribution of the flow field velocity around the floating bodies at the resonance response case of k = 1.55. It is seen that the magnitude of flow speed around Floater 1 is greater than that of Floater 2 in the twin-floater system. The presence of Floater 2 introduces additional fluid disturbances and interference effects in the flow field. These effects superimpose on the flow field originally generated by Floater 1. This results in an increase in the velocity amplitude around Floater 1 in the twin-floater system compared to the case of a single float.

Figure 15: Distribution of the flow velocity field around the floating bodies under resonance condition of k = 1.55.

4.1.2 Force Responses

Fig. 16 shows the variation in the horizontal and vertical force amplitude on the floating bodies without liquid tank under different incident wave cases. It can be observed that the horizontal force on the floating body reaches the peak value at the resonant frequency of wave number k = 1.55. Moreover, a second resonance peak appears at k = 1.9, although this peak is lower than the first peak at k = 1.55. When it comes to the vertical force, the largest peak occurs at k = 1.9, which is much greater than the peak at k = 1.55. The force amplitude of Floater 2 is much lower than that of Floater 1 due to the wave shielding effect.

Figure 16: Force amplitudes of the floating bodies without tank at different wave numbers.

Fig. 17 and Fig. 18 compare the time series of force between the three floaters for the two cases with the wave number k = 1.55 and k = 2.2, respectively. It can be seen that the force signal is not as smooth as the motion signal and the force signal contains high-frequency vibrations. There is a phase difference between the two floaters in the twin-floater system.

Figure 17: Comparison of the force time series of the floating bodies without tank in the case of k = 1.55.

Figure 18: Comparison of the force time series of the floating bodies without tank in the case of k = 1.9.

It can be seen that, compared to the motion signal presented in Fig. 13, the force signal appears less smooth, indicating the presence of high-frequency vibration components in the force signal. When subjected to wave impacts, the pressure distribution on the surface of the floating body is nonlinear. Moreover, there is a sequential order in the occurrence of impact loads along the direction of wave propagation (see Fig. 19). Since force is obtained by integrating pressure over the surface area of the floating body, the force signal exhibits some high-frequency characteristics. Upon further comparison with Floater 1, it is evident that the high-frequency components in the twin-floater system are significantly more numerous than those in the single floater system. Fig. 20 shows the frequency components of the dimensionless force obtained by applying FFT on the time series. Within the primary frequency range of 0–6 Hz, multiple high-frequency fluctuations are present in the horizontal force, whereas only two major high-frequency fluctuations exist in the vertical force. This characteristic makes the high-frequency components of the horizontal force more prominent compared to those of the vertical force. For the high-frequency components of the horizontal force, the energies at 1.1 Hz, 2.1 Hz, and 4 Hz for Floater 1 in the twin-floater system are significantly higher than those in the single floater system. This low frequency 1.1 Hz can be regarded as the main excitation frequency caused by the wave, which is related to the frequency of the wave itself. While the frequency 2.1 Hz is nearly twice the main frequency, it is the frequency-doubling effect of the wave. And the larger 4 Hz causes higher vertical vibration, which may be related to the vertical vibration frequency of the floating body itself. These high-frequency components are the fundamental reason why the time history response of Floater 1 in the twin-floater system is higher than that in the single floater system.

Figure 19: Distribution of the pressure field around the floating bodies under resonance condition of k = 1.55.

Figure 20: Comparison of the frequency component of the floating bodies without tank in the case of k = 1.55. (a) Horizontal force; (b) Vertical force.

4.2 Comparison between One- and Twin-Floater with Liquid Tank

In this part, the flow characteristics and dynamic behavior of the free-floating bodies with internal liquid tank under wave excitation are investigated. The response is compared between the one- and twin-floater systems.

4.2.1 Motion Responses

The variation in the motion amplitude of the floating bodies with liquid tank under different incident wave cases is shown in Fig. 21. The heave amplitude generally exhibits a decrease trend as the wave number increases from k = 1.4 to 2.2. Compared to conventional floating bodies without liquid tanks, those with liquid tanks do not exhibit prominent resonance phenomena. This indicates that the presence of liquid inside the structure can effectively prevent the occurrence of resonance. For the twin-floater system, the motion amplitude of Floater 2 is much lower than that of Floater 1. This trend is consistent with the behavior of twin-floater system without liquid tanks, primarily due to the shielding effect of the overhead structure. In the range of lower wave numbers, the heave amplitude of Floater 1 in the single floater system is lower than that in the twin floater system; however, in the range of higher wave numbers, the heave amplitude of Floater 1 in the single floater system exceeds that in the twin floater system. Nevertheless, in terms of pitch amplitude, Floater 1 in the single and twin systems exhibits opposite patterns. This suggests that the sloshing of the internal liquid tank enhances the heave amplitude response at low frequencies but reduces the pitch response; whereas at high frequencies, the situation is reversed.

Figure 21: Motion amplitudes of the floating bodies with tank at different wave numbers.

Fig. 22 presents the time series of the three floating bodies’ motion in the case of wave number k = 1.7. It is seen that the heave signal shows standard sine in shape, but the pitch signal shows obvious nonlinearity due to internal tank sloshing. The difference in pitch of Floater 1 between single and twin systems is obvious, especially at the troughs. The pitch value of Floater 1 is generally greater than that of Floater 2 due to the obvious trim state of Floater 1 in waves. The motion amplitude of Floater 2 is lower due to the sheltering effect [41].

Figure 22: Comparison of the motion time series of the floating bodies with tank in the case of k = 1.7.

Fig. 23 presents the flow pattern and velocity distribution inside and outside the floating bodies in an resonant state of k = 2.2. It is found that the wave pattern and speed distribution around Floater 1 is the almost the same for the single and twin systems. The flow speed around Floater 2 is much lower than that of Floater 1. And the internal tank sloshing effect in Floater 2 is not significant.

Figure 23: Distribution of the flow velocity field around the floating bodies with tank under condition of k = 2.2.

4.2.2 Force Responses

Fig. 24 shows the variation in the horizontal and vertical force amplitude on the floating bodies with liquid tank under different incident wave cases. Two peaks are identified in the horizontal force curves, which are around k = 1.7 and 1.9. However, only one peak around k = 1.7 can be identified in the vertical force curve, which is different from the results without tank.

Figure 24: Force amplitudes of the floating bodies with tank at different wave numbers.

Fig. 25 illustrates the time history of the dynamic response experienced by the floating bodies with tank under the wave number k = 1.95. Significant differences in phase and amplitude among the three floating bodies can be observed. Further analysis, based on the pressure distribution presented in Fig. 26, reveals that the overall pressure levels of the three floating bodies are similar. However, the free surface shapes within their internal liquid tanks vary. This inconsistency in free surface shapes is a crucial factor contributing to the differences in forces acting on the floating bodies, with the relative positions and morphological changes of the internal and external submerged surfaces playing a decisive role.

Figure 25: Comparison of the force time series of the floating bodies with tank in the case of k = 1.95.

Figure 26: Distribution of the pressure field around the floating bodies under resonance condition of k = 1.95.

To gain a deeper understanding of the frequency components of these dynamic characteristics, Fig. 27 provides an analysis of the frequency features within the time history signals. Within the frequency range of 0 to 6 Hz, the horizontal force components exhibit two distinct high-frequency fluctuation characteristics, while the vertical force components show two significant high-frequency fluctuations only in the case of a single floater with tank. This result indicates that high-frequency components are more prominent in the horizontal force components compared to the vertical force components. For the high-frequency components of the horizontal force, in the twin-floater system, the energy contributed by Floater 1 is significantly higher than that in the single-floater system. This suggests that the fluid dynamic interference between floats in the twin-floater system enhances the energy transfer and accumulation effects at specific frequencies.

Figure 27: Comparison of the frequency component of the floating bodies with tank in the case of k = 1.95. (a) horizontal force; (b) vertical force.

4.3 Effects of the Liquid Tank on the Hydrodynamic Characteristics

Fig. 28 compares the motion time series between the twin floaters with and without internal tank for the case of k = 1.95. It is seen that the heave motion for the floaters with and without internal tank is close although the floater with tank has a lower mean value. However, the pitch of the floater with the tank is much lower than the one without tank. This indicates that the internal tank has the obvious effect on reducing the pitch motion of the floating body, which acts as an anti-rolling tank. A comprehensive comparison of the movement of the floating bodies and fluid field within one wave period between the four schemes is shown in Fig. 29.

Figure 28: Comparison of the motion time series of the floating bodies with and without tank in the case of k = 1.95.

Figure 29: Figure 29: Comparison of movement of the floating bodies and fluid field within one wave period for the case of k = 1.95: (a) Scheme 1, (b) Scheme 2, (c) Scheme 3, and (d) Scheme 4.

Fig. 30 compares the force time series between the twin floaters with and without internal tank for the case of k = 1.95. It is seen that the horizontal force for the floaters with and without internal tank is close. However, the vertical force of the floater with the tank is much lower than the one without the tank. This is due to the fact that the fluid force on the floaters with internal tank contains the resultant forces of both internal and external forces.

Figure 30: Comparison of the force time series of the floating bodies with and without tank in the case of k = 1.95.

5 Conclusions

This paper conducts numerical simulations on the hydrodynamic behaviour of a 2D twin-floater system with and without the liquid tank using the open-source code DualSPHysics. The following conclusions can be obtained:

• (1)The present SPH method well simulates the resonance response of the gap flow between two bodies. The wave elevation height in the narrow gap between two fixed bodies can be five times larger than the incident waves.
• (2)The motion response of the downstream Floater 2 in the twin-floater system is lower than that of the upstream, which is caused by the incident wave shielding effects. Two response peaks can be identified from the amplitude value curves at different wave numbers.
• (3)The internal tank can largely reduce the pitch motion of the floating body while it has little effect on the heave motion. The internal tank can largely reduce the vertical force on the floating body while it has little effect on the horizontal force.

The study in the work is limited to the 2D plane problem of both the fluid field and the floating body; the effect of 3D on the results can be studied in future work. A systematic study on the gap distance on the hydrodynamic behaviour of both multi-floating bodies’ motions and internal sloshing can be further explored. The developed method can be applied to the simulation of real application scenarios, such as the Floating Liquefied Natural Gas (FLNG) system offloading with an LNG ship in a side-by-side arrangement system.