Hydrodynamic Analysis of Twin-Hull Structures Supporting Floating PV Systems in Offshore and Coastal Regions

Abstract

In this paper, a novel model based on the boundary element method (BEM) is presented for the hydrodynamic analysis of floating twin-hull structures carrying photovoltaic panels, supporting the study of wave responses and their effects on power performance in variable bathymetry regions. The analysis is restricted to two spatial dimensions for simplicity. The method is free of any mild-slope assumptions. A boundary integral representation is applied for the near field in the vicinity of the floating body, which involved simple (Rankine) sources, while the far field is modeled using complete (normal-mode) series expansions that are derived using separation of variables in the constant depth half-strips on either side of the middle, non-uniform domain, where the depth exhibited a general variation, overcoming a mild bottom-slope assumption. The numerical solution is obtained by means of a low-order panel method. Numerical results are presented concerning twin-hull floating bodies of simple geometry lying over uniform and sloping seabeds. With the aid of systematic comparisons, the effects of the bottom slope and curvature on the hydrodynamic characteristics (hydrodynamic coefficients and responses) of the floating bodies are illustrated and discussed. Finally, the effects of waves on the floating PV performance are presented, indicating significant variations of the performance index ranging from 0 to 15% depending on the sea state.

1. Introduction

The energy yield of floating photovoltaics (FPVs) is in the spotlight, as offshore photovoltaic (PV) installations present significant advantages over corresponding onshore ones (see [1,2]). These, among others, include the ample surface available for arrangements in farms, the nearshore/coastal regions and the open sea, including locations that are already licensed for offshore wind parks (in the area between wind turbines), as well as the potential of hybridization with offshore wind energy. The development of offshore FPV parks is particularly important for southern European regions, e.g., in the Mediterranean Sea, since solar radiation in southern latitudes is relatively high, nearly 150–200% greater than that of the Atlantic Sea Ocean, the North Sea and Baltic Sea regions [3], while wind and wave potentials are comparatively low. Furthermore, offshore PV installations present increased efficiency due to the cooling effects of water and wind, which are triggered by the interaction of airflow with the solar panels [4]. It is worth mentioning here that, according to a recent report from DNV GL, it is expected that offshore FPVs will reach maturity by 2030 (see also DNVGL-RP-0584-Edition 2021-03).

On the other hand, although several solar farms have been developed on closed water basins, such as lakes, reservoirs and dams, implementing installations in the open sea is a challenging task, as their interaction with several environmental factors is not yet fully understood [5]. In the offshore and nearshore region, safety and viability require the design and construction of resilient FPV structures that can withstand the harsh marine environment. Regarding the deployment of floating structures of relatively large dimensions in nearshore and coastal areas, it is also expected that bathymetric variations will have significant effects on their responses under wave loads, which also affect the performance of the power output due to oscillatory motions of the structure and the panels arranged on the deck. Stability requirements, in conjunction with a lightweight structure with a center of gravity at a relatively increased height above the keel, led to the consideration of a twin hull structure with a more complicated response pattern and resonance characteristics. In this study, a hydrodynamic model is developed to predict the dynamic responses of a floating structure supporting photovoltaic panels on a deck while being subject to wave loads. For the treatment of complicated resonance phenomena, as well as the effects of finite and possibly variable bathymetry, which characterizes nearshore and coastal regions, a general model based on boundary element methods (BEM) is developed, which is capable of modeling the involved phenomena. The model is then systematically applied in selected examples to produce preliminary results, which are illustrative of the effect of dynamic motions on the energy efficiency of a floating unit.

The interaction of free surface gravity waves with floating bodies at intermediate depths in areas that are characterized by non-uniform seabed topographies is a mathematically interesting problem, which can be used to analyze a wide range of applications, such as the design and performance evaluation of ships and other floating structures operating in nearshore areas. Theoretical aspects of the problem of small-amplitude water waves propagating in a finite water depth and their interaction with floating and/or submerged bodies have been presented under various geometric assumptions by many authors [6,7,8] regarding the existence of trapped modes in a channel with obstructions. Furthermore, shallow-water conditions are frequently encountered in marine applications. When floating structures or docks are moored in shallow-water areas, accurate predictions of the motions induced by the prevailing sea state are needed, not only for optimizing the mooring system, depending on the stability needs of each configuration, but also for ensuring that the under-keel clearance remains sufficient for the structure to avoid grounding in extreme (for the area under study) weather and sea conditions. In many applications, the water depth is assumed to be constant, which is practically valid in cases where there are small depth variations or the floating body’s dimensions are small compared to the bottom variation length. However, in applications involving the utilization of floating bodies in coastal waters, the variations of bathymetry may cause significant effects on the hydrodynamic behavior of ships and structures, especially concerning the wave-induced responses. Under the assumption of slowly varying bathymetry, mild-slope model have been developed for the analysis of wave-induced floating body motion [9]. To treat environments that are characterized by steeper bathymetric variations, e.g., near the coast or the entrances of ports and harbors, extended models are required (see, e.g., Ohyama and Tsuchida [10]).

In the present work, a novel method based on BEM is used for the hydrodynamic analysis of twin-hull floating structures with PV systems; the method is capable of treating the effects of varying bathymetry without any mild-slope assumptions. In particular, a low-order panel based on linear wave theory is developed and verified. Following the hybrid formulation by Yeung [11], the present method utilizes the simplicity of Rankine sources, in conjunction with appropriate representations of the wavefield in the exterior semi-infinite domain, as presented by Nestegard and Sclavounos [12] for 2D radiation problems in deep water and by Drimer and Agnon [13] in the case of finite water depth. Τhe far field is modeled using complete (normal-mode) series expansions, which are derived using separation of variables in the two constant-depth half-strips separating the variable bathymetry region from the regions of wave incidence and wave transmission (see Figure 1). Numerical results are presented concerning twin-hull floating bodies of simple geometry over uniform and sloping seabeds. With the aid of systematic comparisons, the effects of bottom slope on the hydrodynamic characteristics (hydrodynamic coefficients and responses) are presented and discussed. Finally, results are presented regarding the effect of wave-induced motions on floating PV performance, indicating significant variations in the performance index ranging from 0 to 15%, depending on the sea state.

2. Mathematical Formulation

The 2D problem concerning the hydrodynamic behavior of a twin hull floating body of arbitrary cross-section in a coastal–marine environment is considered, as illustrated in Figure 1. A Cartesian coordinate system $x=(x_1,x_2,x_3)$ is introduced, with the origin placed at the mean water level, coinciding with the structure’s center of flotation, with the $x_3$-axis pointing upwards. The configuration is considered unchanged in the $x_1$-direction and, therefore, the analysis is limited to the $x_2{x}_3$ plane, modeling a two-dimensional cross-section.

The environment comprises a water layer bounded by the free surface at $x_3=0$ and the rigid bottom at depth $h=h(x_2)$. It is assumed that $h=h(x_2)$ exhibited a general variation, i.e., the corresponding bathymetry is defined by parallel, straight contours lying between two regions of constant but different water depths: $h=h_a$ in the region of wave incidence and $h=h_b$ in the region of transmission. The fluid is assumed to be homogeneous, inviscid and incompressible and its motion irrotational with a small width. The wavefield in the region is excited by a harmonic incident field, with propagation direction normal to the depth contours (along the $x_2$-axis). Without loss of generality, a left-incident wave field is assumed (see Figure 1). Thus, in the context of linearized wave theory, the fluid motion is fully described by the 2D wave potential $\Phi \left(x_2,x_3;t\right)$, with the velocity field being equal to $v\left(x,t\right)=\nabla \Phi \left(x,t\right).$ Assuming that the free-surface elevation and the wave velocities are small, the potential function $\Phi \left(x_2,x_3;t\right)$ satisfies the linearized wave equations (see, e.g., [14,15]). Under these assumptions, the wavefield is time-harmonic and its potential function can be represented by the time-independent (normalized) complex potential function $\phi$ as:

$$\Phi \left(x_2,x_3;t\right)=\mathrm{Re}\left\{-\frac{igA}{\omega }\phi \left(x_2,x_3;\mu \right)\cdot \exp \left(-i\omega \hspace{0.17em}t\right)\right\},$$

(1)

where $H=2A$ is the incident wave height, $g$ is the acceleration of gravity, $\mu ={\omega }^2/g$ is the frequency parameter and $i=\sqrt{-1}$. The free surface elevation is obtained in terms of the wave potential at $x_3=0$ as follows:

$$\eta \left(x_2;t\right)=-\frac{1}{g}\frac{\partial \Phi \left(x_2,0\hspace{0.17em};t\right)}{\partial t}.$$

(2)

In addition to the physical boundaries (floating body, free surface, seabed), we further introduce two vertical interfaces on either side of the body, serving as incidence/radiation/transmission boundaries. Therefore, the boundary $\partial D$ of the two-dimensional domain $D$, occupied by the fluid, is decomposed into eight subsections $\partial D_i,i=1,2,\dots ,8,$ as illustrated in Figure 1, so that $D$ is enclosed by the curve $\partial D={\displaystyle {\cup }_{i=1}^8\partial D_i}$, with $\partial D_1$ and $\partial D_3$ being the right- and left-hand sides, respectively, of the twin-hull’s wetted surface. The sections of $\partial D$ numbered 2, 4 and 8 correspond to the water-free surface, while $\partial D_6$ is the impermeable seabed. Finally, the wave incidence occurs via $\partial D_5$, which also serves, along with $\partial D_7,$ as a radiation interface for the diffracted field due to the presence of the (fixed) body, as well as the radiation fields that develop due to the wave-induced body’s motions.

Apart from the non-uniform domain $D$ containing the floating body, the total flow field is considered to be of infinite length and, therefore, also comprises the uniform semi-infinite subdomains $D_L$ and $D_R$, where the depth is constant and equal to $h_a$ and $h_b$, respectively. Hence, the function $h=h(x_2)$ is of the form:

$$h\left(x_2\right)=\left\{\begin{array}{ll}{h}_a,\hfill & \hfill x_2\le a\\ h\left(x_2\right),& \hfill a<x_2<b\\ h_b,& \hfill x_2\ge b\end{array}\right.$$

(3)

The function $\phi =\phi \left(x_2,x_3;\mu \right)$ appearing in Equation (1) is the normalized potential in the frequency domain, which will hereafter be simply written as $\phi \left(x_2,x_3\right)$. Using standard floating body hydrodynamic theory [15,16], the potential is decomposed as follows:

$$\phi \left(x_2,x_3\right)={\phi }_p\left(x_2,x_3\right)+\frac{\mu }{A}{\displaystyle \sum _{k=2}^4{\xi }_k{\phi }_k\left(x_2,x_3\right)}$$

(4)

where ${\phi }_p\left(x_2,x_3\right)={\phi }_I\left(x_2,x_3\right)+{\phi }_D\left(x_2,x_3\right)$ is the propagating field, with ${\phi }_I\left(x_2,x_3\right)$ being the incident field, which corresponds to the solution of the wave propagation problem across the non-uniform subdomain in the absence of the floating structure and ${\phi }_D\left(x_2,x_3\right)$ being the diffraction potential, which accounts for the presence of the body, fixed in its mean position. Moreover, the functions ${\phi }_k\left(x_2,x_3\right),k=2,3,4$, denote the radiation potentials associated with the motion of the twin-hull structure, corresponding to its three degrees of freedom (DOF), i.e., the linear transverse motion (sway: k = 2), the linear vertical motion (heave: k = 3) and the rotation about the longitudinal ($x_1$) axis (roll: k = 4). Finally, ${\xi }_k,k=2,3,4$, stand for the complex amplitudes of the corresponding wave-induced motions.

The sub-problems, whose solutions provide the potential functions ${\phi }_k\left(x_2,x_3\right),k=p,2,3,4$, in the variable bathymetry region, were formulated as radiation-type problems in the bounded subdomain $D$, with the aid of the following general representations of the wave potential $\phi \left(x_2,x_3\right)$ in the left- and right-side semi-infinite strips $D_L$ and $D_R$, which are obtained using separation of variables (see, e.g., [17]):

$$\begin{array}{l}{\phi }_p^{(L)}\left(x\right)=\left[\exp \left(i{k}_0^{\left(L\right)}{x}_2\right)+C_0^{\left(L\right)}\exp \left(-i{k}_0^{\left(L\right)}{x}_2\right)\right]Z_0^{\left(L\right)}\left(x_3\right)+\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\displaystyle \sum _{n=1}^{\infty }{C}_n^{\left(L\right)}}\exp \left[k_n^{\left(L\right)}\left(x_2-a\right)\right]Z_n^{\left(L\right)}\left(x_3\right),x\in D_L\end{array}$$

(5a)

$$\begin{array}{l}{\phi }_k^{(L)}\left(x\right)=C_0^{\left(L\right)}\exp \left(-i{k}_0^{\left(L\right)}{x}_2\right)Z_0^{\left(L\right)}\left(x_3\right)+\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\displaystyle \sum _{n=1}^{\infty }{C}_n^{\left(L\right)}}\exp \left[k_n^{\left(L\right)}\left(x_2-a\right)\right]Z_n^{\left(L\right)}\left(x_3\right),x\in D_L,k=2,3,4\end{array}$$

(5b)

$$\begin{array}{l}{\phi }_k^{(R)}\left(x\right)=C_0^{\left(R\right)}\exp \left(i{k}_0^{\left(R\right)}{x}_2\right)Z_0^{\left(R\right)}\left(x_3\right)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\displaystyle \sum _{n=1}^{\infty }{C}_n^{\left(R\right)}{Z}_n^{\left(R\right)}}\left(x_3\right)\exp \left[k_n^{\left(R\right)}\left(b-x_2\right)\right],x\in D_R,k=p,2,3,4\end{array}$$

(5c)

The first term $\left(n=0\right)$ in the series Equations (5) is the propagating mode, while the remaining ones $\left(n=1,2,\dots \right)$ are the evanescent modes with $C_n^{(i)}\left(n=0,1,2,\dots \right)$ being the corresponding coefficients. The first term of ${\phi }_p^{(L)}\left(x\right)$ is further separated into a unit-amplitude mode propagating toward D, playing the role of the incident field, and the additional mode $C_0^{\left(L\right)}\exp (-i{k}_0^{\left(L\right)}{x}_2)Z_0^{\left(L\right)}(x_3)$, propagating toward $-\infty$ in the $x_2$-direction, which is the reflected field coming from the diffraction potential ${\phi }_D\left(x\right).$ In the above expansions, the functions ${\left\{Z_n^{(i)}\right\}}_{n=0}^{\infty }$ are defined as $Z_n^{\left(i\right)}=\cosh \left[k_n^{(i)}\left(z+h^{(i)}\right)\right]/\cosh \left(k_n^{(i)}{h}^{(i)}\right)$ and are obtained using separation of variables via the vertical Sturm–Liouville problem, to which Laplace’s equation reduces in the constant depth strips $\left\{\left.D_L\right|-h_a<x_3<0,x_2<a\right\}$ and $\left\{\left.D_R\right|-h_b<x_3<0,x_2>b\right\}$. The corresponding eigenvalues $k_0^{(i)}$ and ${\left\{k_n^{(i)}\right\}}_{n=1}^{\infty }$ are respectively obtained as the real root and the imaginary roots of the dispersion relation: ${\omega }^2=k^{(i)}g\cdot \tanh \left(k^{(i)}{h}^{(i)}\right)\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}i=L,R,$ where $g$ denotes the acceleration due to gravity. The completeness of the expansions derives from the standard theory of regular eigenvalue problems (see, e.g., [18]). Based on the above representations, the hydrodynamic problems concerning the propagating and radiation potentials ${\phi }_k(x_2,x_3)$ were formulated as radiation-type problems, satisfying the following systems of equations, boundary conditions and matching conditions for $\hspace{0.17em}k=p,2,3,4$:

$$\frac{{\partial }^2{\phi }_k(x_2,x_3)}{\partial x_2{}^2}+\frac{{\partial }^2{\phi }_k(x_2,x_3)}{\partial x_3{}^2}=0,x\in {\left.D\right|}_{\left(DomainofTransmission\right)}$$

(6a)

$$\frac{\partial {\phi }_k\left(x\right)}{\partial n}-\mu {\phi }_k\left(x\right)=0,\mu =\frac{{\omega }^2}{g},{\left.x\in \left(\partial D_2\cup \partial D_4\cup \partial D_8\right)\right|}_{\left(FreeSurface\right)}$$

(6b)

$$\frac{\partial {\phi }_k\left(x\right)}{\partial n}=0,x\in {\left.\partial D_6\right|}_{(Seabed)}$$

(6c)

$$\frac{\partial {\phi }_k\left(x\right)}{\partial n}=N_k\left(x\right),{\left.x\in \left(\partial D_1\cup \partial D_3\right)\right|}_{\left(WettedSurface\right)}$$

(6d)

$$\frac{\partial {\phi }_k\left(x\right)}{\partial n}-T_L\left[{\phi }_k^{(L)}\left(x\right)\right]=Q_k,\hspace{0.17em}\hspace{0.17em}x\in {\left.\partial D_5\right|}_{\left(\mathrm{Incidence}/\mathrm{Reflection}/\mathrm{Radiation}\right)}$$

(6e)

$$\frac{\partial {\phi }_k\left(x\right)}{\partial n}-T_R\left[{\phi }_k^{(R)}\left(x\right)\right]=0,\hspace{0.17em}\hspace{0.17em}x\in {\left.\partial D_7\right|}_{\left(\mathrm{Radiation}\right)}$$

(6f)

The above boundary sections are also illustrated in Figure 1. Moreover, in Equations (6a)–(6f), $n=\left(0,n_2,n_3\right)$ denotes the unit normal vector to the boundary $\partial D$, directed to its exterior. The boundary data $N_k,k=2,3,4$ appearing on the right-hand side of Equation (6d) are defined by the components of the generalized normal vector on the wetted surface boundary section $\partial D_1\cup \partial D_3$: $N_2=n_2$, $N_3=n_3$ and $N_4=x_2{n}_3-x_3{n}_2$, and constitute the (unit-amplitude) excitations of the system in Equation (6a)–(6f) for each k. $N_p$ is set to 0 so that the solution of the propagating field is obtained by treating the floating body as an impermeable, immobile solid boundary. Finally, the operators $T_L\left[{\phi }_k^{(L)}\left(x\right)\right]$ and $T_R\left[{\phi }_k^{(R)}\left(x\right)\right]$ are appropriate Dirichlet-to-Neumann (DtN) maps (see, e.g., [19]), ensuring the complete matching of the fields ${\phi }_k\left(x\right),k=p,2,3,4$, on the vertical interfaces $\partial D_5$ and $\partial D_7$, respectively. These operators are derived from Equations (5a)–(5c), exploiting the completeness properties of the vertical bases $\left\{Z_n^{(i)}\left(z\right),n=0,1,2,\dots \right\},i=L,R$. More details are provided in Appendix A.

3. The BEM for Floating Twin-ull Structures

3.1. The Incidence, Diffraction and Radiation Problems

The corresponding problems of the propagating and radiation potentials ${\phi }_k\left(x\right),k=p,2,3,4$, given as Equations (6) were treated by means of boundary integral equation formulations that are based on the single-layer potential (see, e.g., [20]). Accordingly, the following integral representations are introduced for ${\phi }_k\left(x\right),k=p,2,3,4$, in the bounded subdomain $D:$

$${\phi }_k\left(x\right)={\displaystyle \underset{\partial D}{\int }{\sigma }_k\left(x΄\right)G\left(x΄,x\right)d\mathcal{\ell }\left(x΄\right)},x=\left(x_2,x_3\right)\in D,x΄\in \partial D,k=p,2,3,4$$

(7)

where $G\left(x΄,x\right)=\ln \left|x΄-x\right|/2\pi$ is the Green’s function of the Laplace equation in 2D free-space; ${\sigma }_k\left(x΄\right)$ is a source/sink strength distribution, defined on the boundary of the bounded subdomain D for each of the four subproblems; and $d\mathcal{\ell }\left(x΄\right)$ denotes the differential element along the boundary $\partial D$ (see, e.g., [21,22]). Based on the properties of the single-layer distributions, the corresponding normal derivatives of the functions ${\phi }_k\left(x\right),k=p,2,3,4$, on the boundary $\partial D$ are given by (see, e.g., [21]):

$$\frac{\partial {\phi }_k\left(x\right)}{\partial n}=-\frac{{\sigma }_k\left(x\right)}{2}+{\displaystyle \underset{\partial D}{\int }{\sigma }_k\left(x΄\right)\frac{\partial G\left(x΄,x\right)}{\partial n}d\mathcal{\ell }\left(x΄\right)},\hspace{0.17em}\hspace{0.17em}\left(x,x΄\right)\in \partial D.$$

(8)

Using the above in Equations (6b)–(6f), we obtain a system of boundary integral equations, with support on the different sections of $\partial D$ for the determination of the corresponding unknown source distribution ${\sigma }_k\left(x\right),x\in \partial D,k=p,2,3,4$, for each of the potential functions ${\phi }_k\left(x\right),k=p,2,3,4$. The final system read as follows for $k=p,2,3,4$:

$$\begin{array}{l}-\frac{{\sigma }_k\left(x\right)}{2}+{\displaystyle \underset{\partial D}{\int }{\sigma }_k\left(x΄\right)\frac{\partial G\left(x΄,x\right)}{\partial n}d\ell \left(x΄\right)+}\\ \hspace{1em}\hspace{1em}-\mu {\displaystyle \underset{\partial D}{\int }{\sigma }_k\left(x΄\right)G\left(x΄,x\right)d\ell \left(x΄\right)}=0,x\in \left(\partial D_2\cup \partial D_4\cup \partial D_8\right),x΄\in \partial D,\end{array}$$

(9a)

$$-\frac{{\sigma }_k\left(x\right)}{2}+{\displaystyle \underset{\partial D}{\int }{\sigma }_k\left(x΄\right)\frac{\partial G\left(x΄,x\right)}{\partial n}d\ell \left(x΄\right)}=0,x\in \partial D_6,x΄\in \partial D,$$

(9b)

$$-\frac{{\sigma }_k\left(x\right)}{2}+{\displaystyle \underset{\partial D}{\int }{\sigma }_k\left(x΄\right)\frac{\partial G\left(x΄,x\right)}{\partial n}d\ell \left(x΄\right)}=N_k\left(x\right),x\in \left(\partial D_1\cup \partial D_3\right),x΄\in \partial D,$$

(9c)

$$\begin{array}{l}-\frac{{\sigma }_k\left(x\right)}{2}+{\displaystyle \underset{\partial D}{\int }{\sigma }_k\left(x΄\right)\frac{\partial G\left(x΄,x\right)}{\partial n}d\ell \left(x΄\right)}+\\ \hspace{1em}\hspace{1em}-T_L\left[{\displaystyle \underset{\partial D}{\int }{\sigma }_k\left(x΄\right)G\left(x΄,x\right)d\ell \left(x΄\right)}\right]=Q_k,x\in \partial D_5,x΄\in \partial D,\end{array}$$

(9d)

$$\begin{array}{l}-\frac{{\sigma }_k\left(x\right)}{2}+{\displaystyle \underset{\partial D}{\int }{\sigma }_k\left(x΄\right)\frac{\partial G\left(x΄,x\right)}{\partial n}d\ell \left(x΄\right)}+\\ \hspace{1em}\hspace{1em}-T_R\left[{\displaystyle \underset{\partial D}{\int }{\sigma }_k\left(x΄\right)G\left(x΄,x\right)d\ell \left(x΄\right)}\right]=0,x\in \partial D_7,x΄\in \partial D,\end{array}$$

(9e)

From the above systems’ solutions ${\sigma }_k,k=p,2,3,4$, the corresponding potential functions ${\phi }_k(x),k=p,2,3,4$ and all quantities associated with them were calculated using Equations (7) and (8) in the bounded subdomain D. The solutions of the system consisting of Equations (9a)–(9e) are obtained numerically by means of a low-order boundary element method based on simple (Rankine) sources (see also [23]). The geometry of the different sections of $\partial D$ is approximated using linear segments on which the source distribution is taken to be piecewise constant. In this case, the boundary integrals in Equations (9a)–(9e) associated with each element’s contribution can be analytically calculated (see, e.g., [24]) and the systems of boundary integral equations reduce to an equal number of algebraic systems, whose unknowns are the vectors ${\left\{{\sigma }_{k\hspace{0.17em}j}\right\}}_{j=1}^M,k=p,2,3,4$ with $M$ being the number of linear boundary elements used to approximate the geometry of $\partial D$.

3.2. Equations of Motion

The total hydrodynamic loads (forces and moment) on the twin-hull structure consist of the Froude-Krylov loads, which are solely due to the undisturbed incident field ${\phi }_I\left(x\right),$ the diffraction loads caused by the pressure field generated by the presence of the fixed floating body and the radiation loads due to the pressure fields of the wavefields “radiated” by the oscillating body. Based on the calculated propagating potential ${\phi }_p(x)$(consisting of the incident and diffraction potentials), the summation of the Froude-Krylov and the diffraction-induced hydrodynamic forces, as well as the corresponding moment $(F_k,\hspace{0.17em}\hspace{0.17em}k=2,3,4)$, are calculated using surface integration, as follows:

$$F_k=i\omega \rho {\displaystyle \underset{\partial D_1\cup \partial D_3}{\int }{\phi }_P\left(x\right)\cdot N_k\left(x\right)d\ell \left(x\right)},\hspace{0.17em}k=2,3,4,x\in \left(\partial D_1\cup \partial D_3\right)$$

(10)

where $\rho$ denotes the fluid (water) density and $N_k,k=2,3,4$, is the generalized normal vector on the wetted surface (also defined in Section 2). Moreover, from the radiation potentials ${\phi }_k\left(x\right),k=2,3,4$, the hydrodynamic coefficients are calculated using:

$${\omega }^2{A}_{kl}+i\omega B_{kl}=i\omega \rho {\Pi }_{k\hspace{0.17em}l},\hspace{0.17em}\hspace{0.17em}l,k=2,3,4,\mathrm{where}$$

(11a)

$${\Pi }_{kl}={\displaystyle \underset{\partial D_1\cup \partial D_3}{\int }{\phi }_l\left(x\right)N_k\left(x\right)d\ell \left(x\right)},\hspace{0.17em}\hspace{0.17em}l,k=2,3,4,x\in \left(\partial D_1\cup \partial D_3\right)$$

(11b)

In the above expressions, $A_{(3\times 3)}$ is the (symmetric) matrix of the added inertial coefficients, which, for each frequency, corresponded to the proportion of the radiation loads in phase with the structure’s acceleration (in the frequency domain). $B_{(3\times 3)}$ is the corresponding matrix of hydrodynamic damping coefficients, which consists of the part of the radiation loads in phase with the structure’s velocity. Details about the definitions of the hydrodynamic forces and coefficients, as well as the system of equations of motion, can be found in [15] or in ship hydrodynamics textbooks (see, e.g., [16,24]). The latter quantities allow us to formulate and solve the equations of motion of the floating body in the inhomogeneous domain. The general form of the equations of motion in the frequency domain for the 2D twin-hull structure considered is:

$$\left\{-{\omega }^2\left[M+A\left(\omega \right)\right]-i\omega \hspace{0.17em}B\left(\omega \right)+C\right\}\xi =F$$

(12)

where $C$ is the hydrostatic restoring forces and moments acting on the structure.

Due to the symmetry of the body with respect to the vertical axis $x_3$, the component $N_3$ of the generalized normal vector is symmetric, while the components $N_k,k=2,4$, are antisymmetric. Assuming that the seabed profile variations do not significantly alter the radiation potentials ${\phi }_k(x),k=2,3,4$, near the floating structure, the potential function ${\phi }_3(x)$ is also symmetric and the functions ${\phi }_k(x),k=2,4$ are antisymmetric. This fact implies that ${\Pi }_{32}={\Pi }_{34}=0$ and ${\Pi }_{23}={\Pi }_{42}=0$. Therefore, the dynamic equations of motion relating to the oscillations of the body are simplified in the following form, where the heaving motion $({\xi }_3)$ is decoupled from the sway and roll motions $({\xi }_2,{\xi }_4)$ of the twin hull:

$$\left[-{\omega }^2\left(M+A_{22}\right)-i\omega B_{22}\right]{\xi }_2-\left({\omega }^2{A}_{24}+i\omega B_{24}\right){\xi }_4=F_2,$$

(13a)

$$\left[-{\omega }^2\left(M+A_{33}\right)-i\omega B_{33}+2\rho g{B}_{(H)}\right]{\xi }_3=F_3,$$

(13b)

$$\left(-{\omega }^2{A}_{42}-i\omega B_{42}\right){\xi }_2+\left[-{\omega }^2\left(I_{44}+A_{44}\right)-i\omega B_{44}+Mg\cdot GM\right]{\xi }_4=F_4,$$

(13c)

where $B_{(H)}$ is the breadth of each individual hull and $GM$ denotes the metacentric height. The total mass equals $M=\rho \cdot \nabla$, referring to unit length in the transverse direction $(kg/m)$, where $\rho$ denotes the fluid’s density and $\nabla$ is the displacement volume of the structure. Moreover, due to the symmetry of the floating structure, its center of buoyancy $(B)$ is located on the vertical line $x_2=0$ and its $x_3$ coordinate is calculated as the center of area of the submerged volume’s cross-section. The center of gravity $(G)$ is also located at $x_2=0$ due to symmetry of the configuration and its $x_3$ coordinate is considered to be located at the waterplane ($x_3=0$). A total radius of gyration per unit length in the transverse direction of $R_G=\left(B_{(T)}-B_{(H)}\right)/2$ is considered about the longitudinal axis $(x_1)$, where $B_{(T)}$ is the total breadth of the twin-hull structure and, therefore, $I_{44}=M{(R_G)}^2.$ The metacentric radius was evaluated as $BM=I/\nabla$, where I is the second moment of area of the waterplane, calculated using applying Steiner’s theorem as:

$$I=2\left[\left(\frac{B_{(H)}{}^3}{12}\right)+B_{(H)}\cdot {\left(\frac{B_{(T)}-B_{(H)}}{2}\right)}^2\right],$$

(14)

which also refers to the unit length in the transverse direction $(x_1)$. Finally, the metacentric height was calculated as $GM=KB+BM-KG$, where K is a reference point with coordinates $(0,x_3)$. The above equations can also be modified to include other external forces, as e.g., mooring forces or spring terms (see, e.g., Section 3.5 of [25]). The solution of the above system (13) provides us with the complex amplitudes of the corresponding motions of the twin hull: ${\xi }_k,k=2,3,4$. Then, the total wave potential is obtained using Equation (4), from which the hydrodynamic pressure is obtained using Bernoulli’s theorem. The wave loads on the floating structure are calculated using pressure integration on the wetted surface $\partial D_1\cup \partial D_3$.

4. Numerical Results

4.1. Comparison with Other Methods and Verification

The results obtained by the previously described numerical model are here compared to previous research for verification purposes. The results concern a twin-hull floating structure whose individual hulls are cylindrical, with a draft equal to the radius, which results in wetted surfaces whose cross-section shapes are semicircles. Numerical results regarding the above configuration were presented by Ohcusu [26] in 1969 and Rhee [27] in 1982 concerning the amplitude ratio of the radiated fields’ wave height away from the body divided by the amplitude of the forced oscillation that excites the field itself in calm water.

Figure 2 illustrates the aforementioned ratio regarding the heave and sway motions in the case of unit amplitude of the twin-hull structure. The wetted surface of each hull is semicircle of radius R in the $x_2{x}_3$ plane, while each of the two semicircles’ centers are at a distance P from the origin, following the notation of Rhee [27]. Therefore, the two centers are at a distance 2P apart and the configuration is defined so that $2P/R=3$. The results concern the radiation fields that propagate in deep water, which is achieved in the present numerical model by setting the depth as a constant and equal to half the wavelength, for each simulated frequency, as calculated from the dispersion relation for deep-water $(\lambda =2\pi g/{\omega }^2)$.

The domain extends to three wavelengths away from the floating body in both directions and the free surface elevation is evaluated by the discrete BEM model at the last free surface boundary element away from the structure (adjacent to the first boundary element of the radiation boundary). The amplitude ratios of Figure 2 are presented as functions of the non-dimensional frequency parameter ${\omega }^{2}R/g.$

Indicative results are illustrated in Figure 3 concerning wave fields generated by unit-amplitude forced oscillations of the twin hull in sway and heave, with the non-dimensional frequency parameter ${\omega }^{2}R/g$ set to 1. The amplitude ratios are equal to 0.992 and 0.520 for sway and heave, respectively, as also shown in Figure 2.

An identical twin-hull structure was studied by Dabssi et al. [28] in 2008 regarding its hydrodynamic coefficients. Figure 3 illustrates the added mass and damping of the floating structure in heave $({\xi }_3)$. The added mass $(A_{33})$ has been divided by the structure’s mass, while the damping coefficient for heave $(B_{33})$ has been divided by the mass times the angular frequency $\omega$ so that all presented quantities are non-dimensional. It is noted that the displacement in this case does not need to be numerically calculated since it equals the sum of volumes of two half-cylinders of radius R that were considered to extend to unit length in the transverse direction $(x_1)$ and therefore is equal to $\pi R^2$. The results of Figure 4 concern the heaving motion of the twin-hull in a finite water depth $h$, where $h/R=2$. The calculated data sets are presented as functions of the non-dimensional wavenumber $kR$.

4.2. Hydrodynamic Analysis of Floating Twin-Hull Structure in Waves

The effect of sloping seabed environments on the hydrodynamic characteristics of a twin hull is here illustrated by considering the case of a structure of non-dimensional total breadth equal to $B_{(T)}/h=2/3$, with the non-dimensional breadth and draft of each hull set to $B_{(H)}/h=T/h=1/10$, where $h$ denotes the mean water depth of the inhomogeneous domain $D$. The individual hulls that make up the twin-hull layout were modeled via the cross-section of a Wigley hull at $x_1=0$, which is given by the analytical relation:

$$x_2=\mp \frac{B_{(H)}}{2}\cdot \left[1-{\left(\frac{x_3}{T}\right)}^2\right].$$

(15)

The configuration is considered to be located in an inhomogeneous region (see Figure 5 and Figure 6). The center of gravity was selected to coincide with the center of flotation. The center of buoyancy $(B)$, which is calculated as the center of area of the submerged volume’s cross-section, is located at $(x_2=0,x_3=-0.375\cdot T)$ and, thus, the non-dimensional metacentric height of this layout is $GM/h=1.179$.

Numerical results are presented in Figure 6 and Figure 7 concerning the hydrodynamic behavior of this floating structure in constant depth and over two linear shoals characterized by (constant) bottom slopes of 10 and 20%, respectively (see Figure 6a). The shoaling environments were achieved using a linear depth reduction of $2h/3$ and $4h/3,$ respectively, over a depth variation distance of $10B_{(T)}$, with the mean water depth of all three domains of transmission being equal to $h$. The results concerning the homogeneous domain were plotted using solid lines, while the results concerning the inhomogeneous transmission domains with bottom slopes of 10 and 20% were plotted using dashed lines and dotted lines, respectively.

In particular, Figure 6b illustrates the normalized hydrodynamic forces as functions of the non-dimensional wavelength $\lambda /h$ for all three considered domains of transmission, where $\lambda =2\pi /k_0$ is the wavelength corresponding to the mean water depth $h$, as obtained through application of the dispersion relation: ${\omega }^2=k_{0}g\cdot \tanh \hspace{0.17em}(k_{0}h).$ The normalization used for the hydrodynamic forces is ${\tilde{F}}_k=F_k/\rho ghA,\hspace{0.17em}k=2,3,$ where $A$ is the incident wave amplitude.

Figure 6c,d depict the twin hull’s response amplitude operators (RAOs) associated with its two linear motions, i.e., sway $({\xi }_2)$ and heave $({\xi }_3)$. The body’s linear responses were normalized as $RA{O}_k={\tilde{\xi }}_k=\left|{\xi }_k\right|/A,k=2,3$. Finally, in Figure 6e,f corresponding results concerning the hydrodynamic coefficients are presented. The matrix $A_{(3\times 3)}$ of added inertial coefficients and the matrix $B_{(3\times 3)}$ of hydrodynamic damping coefficients were normalized as:

$$\tilde{A}=\frac{A}{\rho }\left(\begin{array}{ccc}{h}^{-2}& h^{-2}& h^{-3}\\ h^{-2}& h^{-2}& h^{-3}\\ h^{-3}& h^{-3}& h^{-4}\end{array}\right),\tilde{B}=\sqrt{\frac{h}{g}}\cdot \frac{B}{\rho }\left(\begin{array}{ccc}{h}^{-2}& h^{-2}& h^{-3}\\ h^{-2}& h^{-2}& h^{-3}\\ h^{-3}& h^{-3}& h^{-4}\end{array}\right).$$

(16)

Figure 7a illustrates the twin hull’s RAO associated with the angular motion, i.e., roll $({\xi }_4)$. The body angular response is normalized as $RA{O}_4={\tilde{\xi }}_4=\left|{\xi }_4\right|/kA,$ with k being the wavenumber corresponding to the mean water depth $h$. The corresponding normalized hydrodynamic moment ${\tilde{F}}_4=F_4/\rho g{h}^{2}A$ is shown in Figure 7b. Figure 7c depicts the diagonal elements $(A_{44},B_{44})$ of the added inertia and damping matrices. Moreover, the (non-diagonal) elements of the symmetric added inertia and hydrodynamic damping matrices are shown in Figure 7d. All results are plotted as functions of the non-dimensional wavelength $\lambda /h$.

Finally, indicative results regarding the total induced wavefields are depicted in Figure 8 for non-dimensional wavelength equal to $\lambda /h=2.4$. In particular, Figure 8a illustrates the real part of the total potential $\phi \left(x\right)$; see Equation (4) for the three considered cases of 0, 10 and 20% bottom slope (see Figure 6a) in the general vicinity of the floating twin-hull structure. Figure 8b depicts the imaginary parts of the corresponding potential functions. The configuration was made dimensional by setting $h=30\hspace{0.17em}\hspace{0.17em}\mathrm{m}$ and an incident field of amplitude $A=1.5\hspace{0.17em}\hspace{0.17em}\mathrm{m}$ has been considered. The alterations to the wavefields due to the non-uniform topography profiles are clearly seen, as the equipotential lines intersect each seabed boundary section perpendicularly, which implies the satisfaction of the impermeability boundary condition.

5. Effects of Floating Structure Response in Waves on Floating PV Performance

The energy efficiency of a floating photovoltaic (FPV) unit is based on several parameters, many of which are the result of the surrounding marine environment. Some of the factors that affect the energy efficiency of FPV are common with corresponding land-based units, with similar power output levels, while others are absent in land installations.

In open seas, there is generally a higher level of humidity than inland, as well as lower ambient temperatures. The decreased temperatures are a result of various factors, which among others, include [29] the water’s transparency, which results in the incoming solar radiation being transmitted to inner layers of the medium rather than just the surface layer, as well as the fraction of incident irradiation that is naturally used for evaporation. Furthermore, the wind speed is usually higher due to long fetch distances compared to land. The above parameters can help to maintain a low operating temperature of the solar cells, which, in turn, leads to close-to-optimal performance of the solar panel. The latter’s efficiency decreases with increasing temperatures. More importantly, PV efficiency is strongly dependent on the angle of incidence (AOI) of solar irradiation, which, in offshore FPV installations, is directly affected by the dynamic wave-induced motions. In particular, the power output of photovoltaic cells is strongly affected by the angle of incidence (AOI) of solar irradiation and the plane of array (POA) irradiance, which is given by the following equation:

$$POA=DNI\cos \left(AOI\right)+DHI+RI,$$

(17)

where DNI, DHI and RI are the direct, diffuse and reflected irradiance components on a tilted surface, respectively. To provide indicative results regarding the effect of wave-induced motions of the floating structure on the power output, we considered an offshore installation in the geographical sea area of the southern Aegean Sea. For the latter area, the optimized values for tilt and azimuth angles of the photovoltaic installations, respectively, are ${\theta }_T={30}^o$ and ${\theta }_A={135}^o$ using data extracted from the Sandia Module Database, which is provided by the PV_LIB toolbox (https://pvpmc.sandia.gov/applications/pv_lib-toolbox/ (accessed on 12 August 2021)).

In this work, a preliminary assessment of a floating photovoltaic system’s energy efficiency is made for twin-hull structures, taking into account data regarding the dynamic motions of the floating unit carrying the panels, as derived by the hydrodynamic model presented earlier, while the interesting effects of temperature and humidity will be studied in future work. The linear motions, i.e., sway (k = 2) and heave (k = 3), are considered to have no important effect on the tilt angle of the panels and, therefore, the angle of incidence. Hence the effect of the unit’s mobility is limited to the angular oscillation i.e., the roll motion (k = 4) under excitation from the beam incident waves.

For this purpose, response data was simulated by assuming specific sea conditions. The latter are characterized by a frequency spectrum used to describe the incident waves. We considered the floating twin-hull structure of total breadth $B_{(T)}=20\hspace{0.17em}\hspace{0.17em}\mathrm{m}$ examined in the previous section in constant water depth $h=30\hspace{0.17em}\hspace{0.17em}\mathrm{m}.$ The sea state is described by a Brettschneider spectrum model (see [30], Chapter 2.3), as follows:

$$S\left(\omega ;H_s,T_p\right)=\frac{5}{16}{H}_s^2\frac{{\omega }_p^4}{{\omega }^5}\exp \left[-\frac{5}{4}{\left(\frac{\omega }{{\omega }_p}\right)}^{-4}\right]$$

(18)

where $H_s$ is the significant wave height, ${\omega }_p=2\pi /T_p$ is the peak frequency and $T_p$ the corresponding peak period.

The roll responses calculated by the present model, as discussed in the previous section, were used to evaluate the fluctuations of the AOI and the effect on the power output performance of a PV system consisting of panels, with the aforementioned values of tilt (relative to the horizontal deck of the structure) and azimuth angles. Specifically, the roll spectrum was calculated using the RAO of the roll motion (see Figure 7a) of the twin-hull structure using:

$$S_4\left(\omega \right)=RA{O}^2\left(\omega \right)k^{2}S\left(\omega \right)$$

(19)

where the wavenumber k is given by the dispersion relation of water waves for the water depth considered. Based on the calculated roll spectrum, time series of the roll motion ${\xi }_4\left(t\hspace{0.17em};H_s,T_p\right)$ of the above floating twin-hull structure were simulated, for the considered configuration (structure and coastal environment) and incident waves, characterized by the parameters $\left(H_s,T_p\right)$ using the random-phase model [30], Chapter 8.2, (see also [31]).

The results were normalized using the value corresponding to calm water (flat horizontal deck of the structure) in the same sea environment, which results in the following definition of the performance index:

$$PI\left(t\right)=\frac{a\hspace{0.17em}\hspace{0.17em}\cos \left({\alpha }_m+{\xi }_4\left(t\right)\right)+b}{a\hspace{0.17em}\hspace{0.17em}\cos \left({\alpha }_m\right)+b}$$

(20)

where $a=DNI$ and $b=DHI+RI$ for the geographical area and sea environment considered, respectively, and ${\alpha }_m$ is a representative value for the angle of incidence.

As an example, the numerical results concerning the calculated roll response of a floating twin-hull structure of breadth $B_{(T)}=20\hspace{0.17em}\hspace{0.17em}\mathrm{m}$ at depth $h=30\hspace{0.17em}\hspace{0.17em}\mathrm{m}$ with an incident wave spectrum (dashed line) and roll angle spectrum (solid line) of the structure in the case of incident waves of significant wave height $H_S=0.5\hspace{0.17em}\hspace{0.17em}\mathrm{m}$ and peak period $T_P=(2\pi /{\omega }_p)=4\hspace{0.17em}\hspace{0.17em}\mathrm{s}$ are presented in Figure 9. Based on the calculated roll spectrum, the simulated time series of the roll motion of the above floating twin-hull structure for the considered coastal environment and incident waves are shown in Figure 10 for a time interval of 1 h. Furthermore, in the same figure, a representative small interval of 3 min was obtained using the random-phase model [30,31] for a sea state that was characterized by significant wave height $H_S=0.5\hspace{0.17em}\hspace{0.17em}\mathrm{m}$ and peak period $T_P=4\hspace{0.17em}\hspace{0.17em}\mathrm{s}$, as generated by winds corresponding to the Beaufort scale levels $BF=1-2$. In this case, indicative results concerning the effect of waves and roll responses of the structure on the performance index are shown in Figure 11, as calculated by Equation (20) using a representative value of the mean angle of incidence ${\alpha }_m=5°$ and omitting, as a first approximation, the effect of diffuse and reflected irradiance components $\left(b\approx 0\right)$. The value of the performance index in calm water was $P{I}_{CALM}=0.9962$. In the considered case of incident waves, which were characterized by a very low energy content, the $RMS$ value of the estimated performance index dropped to $P{I}_{RMS}=0.9947$. The latter’s mean value, as well as the corresponding calm-water value, are shown in Figure 11 using cyan and red lines, respectively.

6. Discussion

Following previous works [32,33], concerning the investigation and modeling of marine renewable energy systems, the present method focused on the estimation of the effect of wave-induced responses on the performance index of a twin-hull FPV structure in various sea conditions, as defined by the wave climatology of the offshore–coastal site where the system was deployed. As an example, the results concerning the performance index that is associated with the wave effects (Equation (20)) of the floating twin-hull structure of breadth $B_{(T)}=20\hspace{0.17em}\hspace{0.17em}\mathrm{m}$ at a water depth $h=30\hspace{0.17em}\hspace{0.17em}\mathrm{m}$ that are presented and discussed above are given in

Table 1. Performance indexes for different sea conditions.

BF	Sea Condition	${\mathit{H}}_{\mathit{S}}\left(\mathit{m}\right)$	${\mathit{T}}_{\mathit{P}}\left(\mathit{s}\right)$	$\mathit{P}{\mathit{I}}_{\mathit{R}\mathit{M}\mathit{S}}$
1–2	1–2	0.5	4	0.9947
3	3	1	6	0.9771
4–5	4	2	8	0.9203
5–6	4–5	3	9	0.8475

for wind waves corresponding to the Beaufort scale from $BF=1$ (relatively calm sea) to $BF=5-6$ conditions.

We observe that the effect of roll responses results in fluctuations of the AOI that could cause a significant drop in the performance index as the sea condition changed from calm to moderate and higher severities. A more complete picture of the sea state’s effect on the FPV module’s power output, as estimated using the present method, is shown in Figure 12, indicating its usefulness for supporting the systematic analysis and design of the system, including the offshore structure, as well as the electric production and storage subsystems. Although inevitable fluctuations of the AOI due to waves in offshore PV units reduce the power output, this negative effect could be balanced or even reversed by the cooling effect and other factors resulting from the marine environment, which is a subject that is left to be considered in future work.

7. Conclusions

A BEM model was developed and applied to the hydrodynamic analysis of twin-hull structures in variable bathymetry regions and was used to predict their hydrodynamic responses and their effects concerning the power output of offshore FPV systems. The analysis was restricted to two spatial dimensions for simplicity. After verification of the method with comparisons against data from the literature, the method was systematically applied and the derived numerical results are presented for floating bodies of simple geometry, lying over uniform and sloping seabeds. With the aid of systematic comparisons, the effects of bottom slope on the hydrodynamic characteristics (hydrodynamic coefficients and responses) of the floating bodies were illustrated and discussed. Finally, response data that was simulated for specific sea conditions, characterized by frequency spectra, were considered to describe the incident waves interacting with a floating twin-hull structure, in order to evaluate the effect of wave-induced fluctuations on the power output performance of the floating PV system. The effects of waves on the floating PV performance are presented, indicating significant variations of the performance index ranging from 0 to 15% depending on the sea state. Future work will be directed toward (a) the detailed analysis of wave and wind environmental factors and their effects on the resulting system’s performance, (b) the extension of the model to 3D including 6-DOF wave motion analysis of the floating structures over general bathymetry and evaluation of their performance and (c) a systematic application of the present method to realistic cases that support the optimized design of floating PV modules in specific marine environments.