**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

**Citation:** Magkouris, A.; Belibassakis, K.; Rusu, E. Hydrodynamic Analysis of Twin-Hull Structures Supporting Floating PV Systems in Offshore and Coastal Regions. *Energies* **2021**, *14*, 5979. https://doi.org/10.3390/en14185979

Academic Editor: Athanasios Kolios

Received: 19 August 2021 Accepted: 16 September 2021 Published: 20 September 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

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–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. The 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.

**Figure 1.** Geometric configuration of the 2D problem.

#### **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*<sup>3</sup> plane, modeling a two-dimensional cross-section.

The environment comprises a water layer bounded by the free surface at *x*<sup>3</sup> = 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* = *ha* in the region of wave incidence and *h* = *hb* 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 leftincident 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 Φ(*x*2, *x*3; *t*), with the velocity field being equal to **v**(**x**, *t*) = ∇Φ(**x**, *t*). Assuming that the free-surface elevation and the wave velocities are small, the potential function Φ(*x*2, *x*3; *t*) 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 *ϕ* as:

$$\Phi(\mathbf{x}\_{2},\mathbf{x}\_{3};t) = \text{Re}\left\{-\frac{i\mathbf{g}\,A}{\omega}\rho(\mathbf{x}\_{2},\mathbf{x}\_{3};\mu)\cdot\exp(-i\omega\,t)\right\},\tag{1}$$

where *H* = 2*A* is the incident wave height, *g* is the acceleration of gravity, *μ* = *ω*2/*g* is the frequency parameter and *<sup>i</sup>* <sup>=</sup> √−1. The free surface elevation is obtained in terms of the wave potential at *x*<sup>3</sup> = 0 as follows:

$$\eta(\mathbf{x}\_2; \mathbf{t}) = -\frac{1}{g} \frac{\partial \Phi(\mathbf{x}\_2, \mathbf{0}; \mathbf{t})}{\partial t}. \tag{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 *∂D* of the two-dimensional domain *D*, occupied by the fluid, is decomposed into eight subsections *∂Di*, *i* = 1, 2, ... , 8, as illustrated in Figure 1, so that *<sup>D</sup>* is enclosed by the curve *<sup>∂</sup><sup>D</sup>* = ∪<sup>8</sup> *<sup>i</sup>*=1*∂Di*, with *∂D*<sup>1</sup> and *∂D*<sup>3</sup> being the right- and left-hand sides, respectively, of the twin-hull's wetted surface. The sections of *∂D* numbered 2, 4 and 8 correspond to the water-free surface, while *∂D*<sup>6</sup> is the impermeable seabed. Finally, the wave incidence occurs via *∂D*5, which also serves, along with *∂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 *DL* and *DR*, where the depth is constant and equal to *ha* and *hb*, respectively. Hence, the function *h* = *h*(*x*2) is of the form:

$$h(\mathbf{x}\_2) = \begin{cases} h\_{a\prime} & \mathbf{x}\_2 \le a \\ h(\mathbf{x}\_2), & a < \mathbf{x}\_2 < b \\ h\_{b\prime} & \mathbf{x}\_2 \ge b \end{cases} \tag{3}$$

The function *ϕ* = *ϕ*(*x*2, *x*3; *μ*) appearing in Equation (1) is the normalized potential in the frequency domain, which will hereafter be simply written as *ϕ*(*x*2, *x*3). Using standard floating body hydrodynamic theory [15,16], the potential is decomposed as follows:

$$
\varphi(\mathbf{x}\_{2}, \mathbf{x}\_{3}) = \varphi\_{p}(\mathbf{x}\_{2}, \mathbf{x}\_{3}) + \frac{\mu}{A} \sum\_{k=2}^{4} \tilde{\varsigma}\_{k} \varphi\_{k}(\mathbf{x}\_{2}, \mathbf{x}\_{3}) \tag{4}
$$

where *ϕp*(*x*2, *x*3) = *ϕI*(*x*2, *x*3) + *ϕD*(*x*2, *x*3) is the propagating field, with *ϕI*(*x*2, *x*3) 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 *ϕD*(*x*2, *x*3) being the diffraction potential, which accounts for the presence of the body, fixed in its mean position. Moreover, the functions *ϕk*(*x*2, *x*3), *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, *ξ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 *ϕk*(*x*2, *x*3), *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 *ϕ*(*x*2, *x*3) in the left- and right-side semi-infinite strips *DL* and *DR*, which are obtained using separation of variables (see, e.g., [17]):

$$\begin{split} \boldsymbol{\Psi}\_{p}^{(L)}(\mathbf{x}) &= \left[ \exp\left(ik\_{0}^{(L)}\mathbf{x}\_{2}\right) + \mathbb{C}\_{0}^{(L)}\exp\left(-ik\_{0}^{(L)}\mathbf{x}\_{2}\right) \right] \mathbf{Z}\_{0}^{(L)}(\mathbf{x}\_{3}) + \\ &\quad + \sum\_{n=1}^{\infty} \mathbb{C}\_{n}^{(L)} \exp\left[k\_{n}^{(L)}(\mathbf{x}\_{2} - a)\right] \mathbf{Z}\_{n}^{(L)}(\mathbf{x}\_{3}), \mathbf{x} \in D\_{L} \end{split} \tag{5a}$$

$$\begin{split} \boldsymbol{\varphi}\_{k}^{(L)}(\mathbf{x}) &= \mathbb{C}\_{0}^{(L)} \exp \left( -i k\_{0}^{(L)} \boldsymbol{x}\_{2} \right) \mathbf{Z}\_{0}^{(L)}(\mathbf{x}\_{3}) + \\ &\quad + \sum\_{n=1}^{\infty} \mathbb{C}\_{n}^{(L)} \exp \left[ k\_{n}^{(L)}(\mathbf{x}\_{2} - a) \right] \mathbf{Z}\_{n}^{(L)}(\mathbf{x}\_{3}), \mathbf{x} \in D\_{L}, k = 2, 3, 4 \end{split} \tag{5b}$$

$$\begin{split} \varphi\_{k}^{(R)}(\mathbf{x}) &= \mathbb{C}\_{0}^{(R)} \exp\left(ik\_{0}^{(R)}\mathbf{x}\_{2}\right) \mathbf{Z}\_{0}^{(R)}(\mathbf{x}\_{3}) \\ &\quad + \sum\_{n=1}^{\infty} \mathbb{C}\_{n}^{(R)} Z\_{n}^{(R)}(\mathbf{x}\_{3}) \exp\left[k\_{n}^{(R)}(b-\mathbf{x}\_{2})\right], \ \mathbf{x} \in D\_{R}, \ k = p,2,3,4 \end{split} \tag{5c}$$

The first term (*n* = 0) in the series Equations (5) is the propagating mode, while the remaining ones (*n* = 1, 2, . . .) are the evanescent modes with *C*(*i*) *<sup>n</sup>* (*n* = 0, 1, 2, . . .) being the corresponding coefficients. The first term of *ϕ*(*L*) *<sup>p</sup>* (**x**) is further separated into a unitamplitude mode propagating toward *D*, playing the role of the incident field, and the additional mode *C*(*L*) <sup>0</sup> exp(−*ik*(*L*) <sup>0</sup> *<sup>x</sup>*2)*Z*(*L*) <sup>0</sup> (*x*3), propagating toward −∞ in the *x*2-direction, which is the reflected field coming from the diffraction potential *ϕD*(**x**). In the above expansions, the functions *Z*(*i*) *n* <sup>∞</sup> *n*=0 are defined as *Z*(*i*) *<sup>n</sup>* <sup>=</sup> cosh *k* (*i*) *n* - *z* + *h*(*i*) / cosh- *k* (*i*) *<sup>n</sup> h*(*i*) and are obtained using separation of variables via the vertical Sturm–Liouville problem, to which Laplace's equation reduces in the constant depth strips { *DL*| − *ha* < *x*<sup>3</sup> < 0, *x*<sup>2</sup> < *a*} and { *DR*| − *hb* < *x*<sup>3</sup> < 0, *x*<sup>2</sup> > *b*}. The corresponding eigenvalues *k* (*i*) <sup>0</sup> and *k* (*i*) *n* <sup>∞</sup> *n*=1 are respectively obtained as the real root and the imaginary roots of the dispersion relation: *<sup>ω</sup>*<sup>2</sup> <sup>=</sup> *<sup>k</sup>*(*i*)*<sup>g</sup>* · tanh- *k*(*i*)*h*(*i*) , *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 *ϕk*(*x*2, *x*3) were formulated as radiation-type problems, satisfying the following systems of equations, boundary conditions and matching conditions for *k* = *p*, 2, 3, 4:

$$\frac{\partial^2 \varrho\_k(\mathbf{x}\_2, \mathbf{x}\_3)}{\partial \mathbf{x}\_2^2} + \frac{\partial^2 \varrho\_k(\mathbf{x}\_2, \mathbf{x}\_3)}{\partial \mathbf{x}\_3^2} = 0, \ \mathbf{x} \in D|\_{\left(Domain\ of\ transmission\right)}\tag{6a}$$

$$\frac{\partial \varphi\_k(\mathbf{x})}{\partial \mathbf{n}} - \mu \varphi\_k(\mathbf{x}) = 0, \; \mu = \frac{\omega^2}{\mathcal{S}}, \; \mathbf{x} \in \left(\partial D\_2 \cup \partial D\_4 \cup \partial D\_8\right)|\_{\text{(Free Surface)}}\tag{6b}$$

$$\frac{\partial \rho\_k(\mathbf{x})}{\partial n} = 0, \; \mathbf{x} \in \left. \partial D\_6 \right|\_{\left( \text{Scaled} \right)} \tag{6c}$$

$$\frac{\partial \rho\_k(\mathbf{x})}{\partial n} = N\_k(\mathbf{x})\_\prime \mathbf{x} \in (\partial D\_1 \cup \partial D\_3)|\_{\text{(Wetted Surface)}}\tag{6d}$$

$$\frac{\partial \rho\_k(\mathbf{x})}{\partial \mathbf{n}} - T\_L \left[ \rho\_k^{(L)}(\mathbf{x}) \right] = Q\_{k'} \mathbf{x} \in \left. \partial D\_5 \right|\_{\text{(Incidence /Reflection /Radiation)}} \tag{6e}$$

$$\frac{\partial \varrho\_k(\mathbf{x})}{\partial n} - T\_{\mathbb{R}} \left[ \varrho\_k^{(\mathbb{R})}(\mathbf{x}) \right] = 0, \; \mathbf{x} \in \left. \partial D\_7 \right|\_{\left( \text{Radiation} \right)} \tag{6f}$$

The above boundary sections are also illustrated in Figure 1. Moreover, in Equations (6a)–(6f), *n* = (0, *n*2, *n*3) denotes the unit normal vector to the boundary *∂D*, directed to its exterior. The boundary data *Nk*, *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 *∂D*<sup>1</sup> ∪ *∂D*3: *N*<sup>2</sup> = *n*2, *N*<sup>3</sup> = *n*<sup>3</sup> and *N*<sup>4</sup> = *x*2*n*<sup>3</sup> − *x*3*n*2, and constitute the (unit-amplitude) excitations of the system in Equation (6a)–(6f) for each *k*. *Np* 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 *TL ϕ*(*L*) *<sup>k</sup>* (**x**) and *TR ϕ*(*R*) *<sup>k</sup>* (**x**) are appropriate Dirichlet-to-Neumann (DtN) maps (see, e.g., [19]), ensuring the complete matching of the fields *ϕk*(**x**), *k* = *p*, 2, 3, 4, on the vertical interfaces *∂D*<sup>5</sup> and *∂D*7, respectively. These operators are derived from Equations (5a)–(5c), exploiting the completeness properties of the vertical bases *Z*(*i*) *<sup>n</sup>* (*z*), *<sup>n</sup>* <sup>=</sup> 0, 1, 2, . . . , *i* = *L*, *R*. More
