Figure 1 illustrates a comprehensive floating energy system, featuring a floating semi-submersible platform with a concave wall, where the top surface can be mounted with solar photovoltaic panels, in addition to four spherical oscillating buoys around it. Each heaving buoy is linked to the central platform through two spring components, allowing for the utilization of suitable power take-off mechanisms to extract energy and convert it into electricity. The floating platform and heaving buoys are motivated by ocean waves, with their mutual dependence realized through the springs connecting them. In this study, we focus solely on analysing the heaving motions of the platform and four buoys, as well as their coupling effects through the spring components.

To better describe the integrated system, a Cartesian coordinate system is established with the origin O at the mass centre of the platform, also coinciding with the undisturbed water surface. The x-axis is along the direction of wave propagation, while the y-axis is perpendicular to it. The z-axis points upwards. The semi-submersible platform is featured as a solid of revolution with a concave wall. The illustration of the generatrix for the central platform is displayed in Figure 2. It has a radius R P 1 at the top and the bottom, and at the waterline, the radius is R P 2 . The curved part is a minor arc of a circle. This design approach enables the achievement of a compact waterplane, while also ensures sufficient space for the placement of components such as springs and buoys. In addition, the mean draught, equalling to the half-height of the platform, is denoted as D.

Four spheroid buoys with same geometry, denoted as B i   i = 1 , 2 , 3 , 4 , are positioned with their average centres evenly spaced on a circle with radius R P 1 on the water surface, or x , y = R P 1 , 0 , − R P 1 , 0 , 0 , − R P 1 and 0 , R P 1 , respectively. The equation of these spheroid is given as x 2 + y 2 / R L 2 + z 2 / R H 2 = 1 , where R L and R H are the equatorial radius of the spheroid and the distance from the centre to the pole along the symmetry axis, respectively. The upper and lower spring components connected to B i are denoted as S u , i and S l , i , respectively, where the subscripts u and l refer to upper and lower. The seabed is assumed to be flat, and the average water depth is denoted as H .

The motions of the heaving buoys B i   i = 1 , 2 , 3 , 4 and the central platform are coupled through the forces due to the spring components, denoted as F s , i i = 1 , 2 , 3 , 4 , and hydrodynamic forces. For each F s , i , it includes the forces due to the upper spring S u , i and that of the lower spring S l , i , which can be written as F s , i = − b u , i + b l , i h ˙ i − h ˙ p − k u , i + k l , i h i − h p (1) where b u , i are the damping coefficients of the upper springs, and b l , i are those of the lower springs, and k u , i and k l , i are the stiffnesses of the upper and lower springs respectively. h i and h p are the vertical displacement of buoy B i and the central platform respectively, the dot above h i and h p denotes the time derivative. Considering the effect of spring systems, the motion equation for each buoy B i can be written as m i h ¨ i + b u , i + b l , i h ˙ i − h ˙ p + k u , i + k l , i h i − h p = F i , i = 1 ∼ 4 (2) where m i is the mass of the buoy B i , and F i is the hydrodynamic force exerted to the buoy B i . The motion of the central platform can be written as m p h ¨ p + ∑ i = 1 4 b u , i + b l , i h ˙ p − h ˙ i + k u , i + k l , i h p − h i = F p , (3) where m p is the mass of the platform, and F p is the fluid force exerted to the platform. Eqs (2) and (3) provide five equations in total, and they should be solved together to obtain the motions of four buoys and the central platform. In the coupling system, the average energy absorbed by one spring component may be expressed as E i = b u , i + b l , i m T × ∫ t t + m T h ˙ i − h ˙ P 2   d t , (4) where m is positive integer and T refers to the period.

We further assume that the motion becomes sinusoidal with time with frequency ω , we may express h i − h p = Re β i e i ω t , where i = − 1 and β i is complex. From (4), we have E i = 1 2 b u , i + b l , i ω 2 ξ i − p 2 . (5)

Here, ξ i − p = β i is the average amplitude of the difference between the instantaneous motions of buoys and platform when the motion becomes stable. Here, we define ξ i and ξ p respectively as the motion amplitude of buoy B i and the platform, we should note that ξ i − p is unequal to ξ i − ξ p due to phase difference of buoys and platform.

Due to the coupling of the motion of buoys and the platform, directly determining the optimal damping coefficients of the springs is challenging. However, these coefficients are crucial for identifying the most effective operational state of the comprehensive system. To address this, we simplify the scenario by using a system comprising one buoy and a spring pair to estimate the optimal damping coefficients for the spring components. In this simplified set-up, the upper and lower ends of the springs, initially connected to the platform, are assumed to be fixed in space. Meanwhile, we allow the single buoy with two springs to move in response to waves. This leads to h p = 0 in Eq. (2), and the motion equation of a single buoy becomes m h ¨ + b u + b l h ˙ + k u + k l h ˙ = F , (6) where h refers to the vertical displacement of a single buoy, and F is the external force from the fluid. Based on linearized velocity potential flow theory, F in (6) comprises wave exciting force f , wave radiation force and restoring force. Removing the radiation and restoring force terms to the left-hand side of Eq. 6, we have m + m a h ¨ + b a + b u + b l h ˙ + k a + k u + k l h = f , (7) in which m a and b a are added mass coefficient and damping coefficient due to wave radiation, respectively. The parameter k a = ρ w g π R P 2 2 refers to the coefficient of restoring force of the buoy, in which g is the acceleration due to gravity. Similarly, we rewrite h = Re β e i ω t and f = Re f 0 e i ω t , where β and f 0 are the complex amplitudes of motion of the buoy and wave exciting force respectively. Substituting these into Eq. 7, the motion amplitude ξ of the single buoy can be obtained as ξ = β = f 0 m + m a ω 2 − k a + k u + k l 2 + b a + b u + b l 2 ω 2 . (8)

From (8), the natural frequency ω n can be obtained as ω n = k a + k u + k l m + m a . (9)

Therefore, the mechanical energy E absorbed by the buoy per second takes the form of E = b f 0 2 2 m + m a ω − k a + k u + k l ω 2 + b a + b 2 , (10) where b = b u + b l . Applying ∂ E ∂ b = 0 , the optimal damping coefficient b opt can be obtained as b opt = m + m a ω 2 − k a + k u + k l 2 ω 2 + b a 2 . (11)

To select the appropriate geometric shape for the buoy, we may utilize Eqs. 10 and (11) for an approximate estimation. In these equations, the hydrodynamic coefficients such as the added mass m a , damping coefficient b a , and the amplitude of wave exciting force f 0 can be obtained through linear analysis using software ANSYS-AQWA. Three spheroids of different parameters are chosen and displayed in Table 1, respectively correspond to an oblate spheroid, a sphere, and a prolate spheroid.

They have the same mass, which is m = 99.129 kg . For simplification, the properties of the spring systems are assumed the same. The damping coefficients of the upper and lower springs are set as b u = b l = 0.5 b opt , where b opt can be calculated from Eq. 11. The stiffnesses of the upper and lower springs are k u = k l = 500 N / m . For these three spheroids, the motion amplitude per wave amplitude and absorbed energy capacity, defined as energy absorbed per second per body mass per wave amplitude, are respectively plotted against wave frequencies in Figures 3, 4, 5 for an oblate spheroid, a sphere and a prolate spheroid. In these three cases, the added mass m a are calculated through linear hydrodynamic analysis and substituted into Eq. 9 to work out the natural frequency ω n , namely, ω n ≈ 6.03  rad / s , 6.05  rad / s and 5.48  rad / s respectively for the oblate spheroid, sphere and prolate spheroid, as marked by the red dashed vertical lines in Figures 3–5.

We may investigate these curves further. According to Eq. 8, it can be observed that if the term b a + b u + b l 2 ω 2 , which relates to the damping force of the spring, is negligible, the motion amplitude ξ would reach its maximum when the wave frequency is equal to the natural frequency, indicating resonance. This is evident in the sphere and prolate spheroid cases, as shown in Figure 4A and 5A, where the greatest motion amplitude is observed near ω = ω n . In addition, as shown in Figure 5B, the maximum energy capacity is also concentrated at the resonance frequency for the prolate spheroid case. However, in the case of the sphere, as illustrated in Figure 4B, there is only a minor energy peak near ω = ω n , while the significant peak appears at ω ≈ 2.5  rad / s . This can also be elucidated by Eq. 10, where if the term [ m + m a ω − k a + k u + k l ω ] 2 and the hydrodynamic damping coefficient b a are both small, then the energy E primarily depends on f 0 2 / b opt . For the sphere, the peak of f 0 2 / b opt occurs at ω ≈ 2.8  rad / s , thus the largest energy capacity is observed near this frequency.

For the oblate spheroid case, as the damping force of spring is noticeably important than other force components, the energy peak mainly depends on the peak of f 0 2 / b opt , the largest energy is at ω ≈ 2.5  rad / s , and the effect of resonance is not prominent in this case, as illustrated in Figure 3. It may be worth noting that, in practice, high absorbed energy in lower wave frequency is expected, and this is for a more energy would be absorbed from the wave with a larger wavelength relative to the size of body. Thus, in the following work, we use the oblate spheroid, with the coefficients R H = 0.228 m and R L = 0.456 m , to consider the coupling effect of the combined system.

The time-domain simulation of wave interaction with the integrated system is conducted. Four oblate spheroid buoys with identical geometries, as outlined in Table 1, along with a floating platform with dimensions R P 1 = 3 m , R P 2 = 0.95 m and D = 2.5 m . The radius of the concave wall of the platform is set as 2.55 m. A cubic computational domain is adopted for the simulation. For each wave frequency, the length of the computational domain is set at 8 λ , and the length of the wave-absorbing zone at both ends of the domain is configured to be 2 λ . In addition, the domain width is fixed at 16m, and the water depth at 10 m.

The convergence study is conducted based on a case study where the circular radian frequency ω of the incident wave is set at 2.512 rad / s , and the corresponding wavelength, denoted as λ , is 9.8 m. The wave height is fixed at H W = 0.15 m . The damping coefficients of the springs are set as b u = b l = b opt / 2 . In the present case, b opt = 2169 N ⋅ s / m . The stiffnesses of the upper and lower springs are set to k u = k l = 500 N / m , as in Section 4.1.

Trimmed meshes are used to discretize the entire computational domain. To account for the motion of bodies, an overset mesh approach is employed. Within the overset domain, the finer mesh relative to the body size is utilized. To enhance numerical efficiency, a multi-block grid has been used. Finer meshes are adopted to regions within the overset domain and in close proximity to the free surface. For example, the region close to the free surface is encrypted four times in z direction. The remaining domain is divided into four subdomains, as illustrated in Figure 6. Mesh sizes within subdomain 1, 2, 3 and 4 are established as 2 l 0 , 4 l 0 , 8 l 0 , and 16 l 0 , respectively. For mesh convergence analysis, two sets of grid lengths are considered as basis, namely, l 0 = 0.05 m and 0.1m, respectively.

In addition, time step (dt) is set at 0.005s and 0.0025s, respectively. Figure 7 presents the studies of the mesh convergence and time step convergence for the vertical motion of platform in the integrated system, while Figure 8 depicts those for buoy 1. Similar convergence results are observed for Buoys 2 to 4 and hence only the results of buoy 1 are provided. The good agreement seen in Figure 7 and 8 verifies that the present procedure is dependent of both time and mesh. Subsequently, simulations were conducted using l 0 = 0.01 m and d t = 0.005 s in the following sections.

In this section, we explore the impact of incident wave frequencies. The parameters of the buoys and the central platform remain consistent with those utilized in the convergency study, despite being subjected to varying wave frequencies. Specifically, six wave frequencies are selected as 0.628 rad / s − 1 , 1.256 rad / s − 1 , 1.884 rad / s − 1 , 2.512 rad / s − 1 , 3.14 rad / s − 1 and 3.768 rad / s − 1 , with a fixed wave height of H W = 0.15 m .

The optimal damping coefficient, b o p t , varying with the wave frequency ω , is calculated using Eq. 9. In the present numerical simulation, the setting of b o p t can be achieved by setting the damping coefficients of the two spring pairs corresponding to B i as b u , i = b l , i = b o p t / 2 . The stiffness of the springs should be set as positive constants, reflecting the capacity to restore energy. In Eq. 13, it is evident that increasing the stiffness raises the natural frequency of the system, which does not favour the absorption of energy in the low frequency region. Therefore, we assign a smaller value to the stiffness, setting k u , i = k p , i = 500 N / m .

Figure 9A illustrates a comparison of the prediction of the motion amplitudes between a single platform analysed in the frequency domain using linearized potential flow theory and the platform coupled with springs and buoys in the integrated system, as analysed in the time domain simulation based on Reynolds-Averaged Navier-Stokes (RANS) equations. The latter one is obtained through a long-time simulation to achieve its steady state. It is observed that the motion amplitudes of platform in the integrated system are generally smaller than those of a single platform in the frequency-domain simulation. This discrepancy may be attributed to several factors. First, the inclusion of buoys alters the mass of the entire system, along with added mass and wave-exciting forces. In addition, a portion of the wave energy is absorbed by the wave energy converters surrounding the platform, leading to a reduction in platform motion. As the wave frequency increases, the motion amplitudes for both cases gradually diminish. It is expected that as the wavelength becomes very short, numerous wave elements will surround the bodies, potentially resulting in the cancellation of the hydrodynamic effects.

In addition to the central platform, Figure 9B illustrates the motion amplitudes of the four buoys in the coupled system, along with the motion amplitude of the single buoy in frequency domain based on linearized analysis. As the symmetry of the positions of B 3 and B 4 corresponding to the incident waves and the central platform, their results coincide with each other. From the figure, we can observe that for wave frequencies other than ω = 2.512  rad / s and 3.14 rad / s , the motion amplitude of B 1 , situated at the upstream side of the platform, is the largest. Those at the two sides are smaller but still larger than that in the frequency-domain prediction, and the one at the downstream. However, at ω = 2.512  rad / s and 3.14 rad / s , the motion amplitude of the single buoy in the frequency-domain prediction is the largest. The motion of the buoy B 1 is smaller than the buoys at the two sides, and at ω = 3.14  rad / s , it is even smaller than that of B 2 at the downstream side.

This observation appears to contradict the intuition that the motion amplitude of buoy B 1 should be largest, as it directly faces the incoming waves, thus experiencing the direct impact of wave energy. While this holds true in many cases, it is important to notice that waves can also be reflected by the central platform. If the incident wave interacts with the reflected wave, the wave energy or wave height in front of the central platform may diminish. Consequently, a significant portion of energy may either propagate to the back through bypassing the platform or forming diffracted waves. Thus, we could see the motion of B 1 is smallest at ω = 3.14  rad / s . Below ω = 2.512  rad / s , the motion amplitudes of buoys B 1 , B 3 and B 4 are larger than those of frequency domain, the reason for this is because among these frequencies, the reflected wave and incident wave are not fully cancelled by each other, they both take positive effects on the motion of buoys.

From the previous analysis and Eq. 5, it becomes evident that the relative motion between the platform and buoy B i , or ξ i − p ( i = 1 , 2 , 3 , 4 ), is crucial for energy absorption. Consequently, Figure 10 illustrates ξ i − p and the corresponding energy capacities, defined as energy absorbed per body mass per wave amplitude, for each buoy B i . In the four wave frequencies other than ω = 2.512  rad / s and 3.14  rad / s , the relative motion of B 1 corresponding to the central platform is larger than that of the other buoys. It remains to be smaller than the motion of the single uncoupled buoy in the frequency-domain prediction. It might be worth noting that ξ 3 exhibits significant differences compared to ξ 2 in these wave frequencies, while ξ 3 − p appears close to ξ 2 − p as shown in Figure 10A. This disparity can be attributed to the phase difference among the buoys.

At ω = 2.512  rad / s and 3.14  rad / s , the relative motion of B 1 , or ξ 1 − p , is no longer the largest. This can be due to the offset between the incident wave and the reflected wave. Overall, the differences in relative motions among buoys are much smaller than those in absolute motions. In other words, while each buoy in the integrated system shows distinct motion characteristics, their motions relative to the central platform are closer due to the coupling effect of the system. The absorbed energy, directly related to the relative motion, follows a similar variation trend. Both relative amplitude and absorbed energy are smaller than the results of a single uncoupled buoy in the frequency domain. This reduction in absorbed energy can be attributed to the change of wave energy distribution caused by the platform within the coupled system.

Figure 11 presents the results obtained with different damping coefficients. Specifically, b p t o is varied as 0, b o p t / 3 , 2 b o p t / 3 and b o p t , which can also be achieved by simply setting b u , i = b l , i = b p t o / 2 . In the present study case, we maintain k u , i = k p , i = 500 N / m , with a wave frequency of ω = 1.256  rad / s and a wave height of H W = 0.15 m .

Figure 11A illustrates the motion of the platform in the frequency-domain analysis without buoys and the time-domain analysis integrated system coupled with buoys, with the latter exhibiting smaller motion response compared to the former. One reason for this discrepancy is the absorption of wave energy by the surrounding buoys, while another factor is the alteration of the system’s dynamic performance due to the increased mass, added mass etc. However, it is noteworthy that the variation in platform motion with changes in the springs’ damping coefficient is not significant. Therefore, we can conclude that the variation in the dynamic performance of the system exerts more profound effects on the motion of the platform. Figure 11B displays the absolute motion amplitudes of the buoys, and the frequency-domain prediction of a single buoy is also displayed for comparison. Among the buoys, B 1 exhibits the largest motion, followed by B 3 and B 4 , and B 2 with the smallest motion. The frequency-domain prediction results for the single buoy fall between B 2 and B 3 or B 4 . These variations can be due to the interaction between incident waves, radiation waves, and diffracted wave, as previously discussed.

The relative motion results are illustrated in Figure 12A, with B 2 , B 3 and B 4 closely aligned, while B 1 exhibits greater movement. This trend is mirrored in the energy distribution displayed in Figure 12B. As the value of b p t o increases, both absolute motion and relative motion decrease, while the absorbed energy increases. At b p t o = b o p t , the energy reaches the peak. This suggests that the optimal damping coefficient obtained from the single buoy based on the frequency-domain analysis can be used to the prediction of the optimal damping coefficient for the integrated system.