**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 *ϕk*(**x**), *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 *ϕk*(**x**), *k* = *p*, 2, 3, 4, in the bounded subdomain *D*:

$$\varphi\_k(\mathbf{x}) = \int\_{\partial D} \sigma\_k(\mathbf{x'}) G(\mathbf{x'}, \mathbf{x}) d\uparrow\_\star(\mathbf{x'}),\\\mathbf{x} = (\mathbf{x}\_2, \mathbf{x}\_3) \in D,\\\mathbf{x'} \in \partial D, \ k = p, 2, 3, 4 \tag{7}$$

where *G*(**x** , **x**) = ln|**x** − **x**|/2*π* is the Green's function of the Laplace equation in 2D freespace; *σk*(**x** ) is a source/sink strength distribution, defined on the boundary of the bounded subdomain *D* for each of the four subproblems; and *d*(**x** ) denotes the differential element along the boundary *∂D* (see, e.g., [21,22]). Based on the properties of the single-layer distributions, the corresponding normal derivatives of the functions *ϕk*(**x**), *k* = *p*, 2, 3, 4, on the boundary *∂D* are given by (see, e.g., [21]):

$$\frac{\partial \rho\_k(\mathbf{x})}{\partial n} = -\frac{\sigma\_k(\mathbf{x})}{2} + \int\_{\partial D} \sigma\_k(\mathbf{x}') \frac{\partial G(\mathbf{x}', \mathbf{x})}{\partial n} d\upzeta(\mathbf{x}'), \text{ (x, x')} \in \partial D. \tag{8}$$

Using the above in Equations (6b)–(6f), we obtain a system of boundary integral equations, with support on the different sections of *∂D* for the determination of the corresponding unknown source distribution *σk*(**x**), **x** ∈ *∂D*, *k* = *p*, 2, 3, 4, for each of the potential functions *ϕk*(**x**), *k* = *p*, 2, 3, 4. The final system read as follows for *k* = *p*, 2, 3, 4:

$$\begin{array}{ll} -\frac{\sigma\_{k}(\mathbf{x})}{2} + \int \sigma\_{k}(\mathbf{x'}) \frac{\partial \mathbb{G}(\mathbf{x'}, \mathbf{x})}{\partial n} d\ell(\mathbf{x'}) + \\ -\mu \int \sigma\_{k}(\mathbf{x'}) \mathbb{G}(\mathbf{x'}, \mathbf{x}) d\ell(\mathbf{x'}) = 0, \ \mathbf{x} \in (\partial D\_{2} \cup \partial D\_{4} \cup \partial D\_{8}), \ \mathbf{x'} \in \partial D, \end{array} \tag{9a}$$

$$-\frac{\sigma\_k(\mathbf{x})}{2} + \int\_{\partial D} \sigma\_k(\mathbf{x'}) \frac{\partial \hat{G}(\mathbf{x'}, \mathbf{x})}{\partial n} d\ell(\mathbf{x'}) = 0, \ \mathbf{x} \in \partial D\_{\delta\prime} \ \mathbf{x'} \in \partial D\_{\prime} \tag{9b}$$

$$-\frac{\sigma\_k(\mathbf{x})}{2} + \int\_{\partial D} \sigma\_k(\mathbf{x'}) \frac{\partial G(\mathbf{x'}, \mathbf{x})}{\partial n} d\ell(\mathbf{x'}) = \mathcal{N}\_k(\mathbf{x}), \ \mathbf{x} \in (\partial D\_1 \cup \partial D\_3), \ \mathbf{x'} \in \partial D\_\prime \tag{9c}$$

$$\begin{split} -\frac{\sigma\_{k}(\mathbf{x})}{2} &+ \int\_{\partial D} \sigma\_{k}(\mathbf{x'}) \frac{\partial G(\mathbf{x'}, \mathbf{x})}{\partial n} d\ell(\mathbf{x'}) + \\ &- T\_{L} \left[ \int\_{\partial D} \sigma\_{k}(\mathbf{x'}) G(\mathbf{x'}, \mathbf{x}) d\ell(\mathbf{x'}) \right] = Q\_{k'} \ \mathbf{x} \in \partial D\_{5}, \ \mathbf{x'} \in \partial D\_{\prime} \\ &- \frac{\sigma\_{k}(\mathbf{x})}{2} + \int\_{\partial D\_{\prime}} \sigma\_{k}(\mathbf{x'}) \frac{\partial G(\mathbf{x'}, \mathbf{x})}{\partial n} d\ell(\mathbf{x'}) + \\ &- T\_{R} \left[ \int\_{\partial D} \sigma\_{k}(\mathbf{x'}) G(\mathbf{x'}, \mathbf{x}) d\ell(\mathbf{x'}) \right] = 0, \ \mathbf{x} \in \partial D\_{7}, \ \mathbf{x'} \in \partial D\_{\prime} \end{split} \tag{9e}$$

From the above systems' solutions *σk*, *k* = *p*, 2, 3, 4, the corresponding potential functions *ϕ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 *∂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 *<sup>σ</sup>k j<sup>M</sup> j*=1 , *k* = *p*, 2, 3, 4 with *M* being the number of linear boundary elements used to approximate the geometry of *∂D*.
