A Novel BEM for the Hydrodynamic Analysis of Twin-Hull Vessels with Application to Solar Ships

The coordinate system $x = (x_{1}, x_{2}, x_{3})$ is used, with the origin placed at Mean Water Level $(MWL)$ , coinciding with the center of flotation. The $x_{1}$ -axis runs parallel to the vessel’s length, $x_{2}$ is the transverse axis, and the vertical coordinate $x_{3}$ is negative in the water body. Let $Ω \in ℝ^{3}$ denote a flow domain that extends infinitely to all azimuthal directions $θ = a t a n 2 (x_{1}, x_{2}),$ as well as the negative $x_{3}$ -direction, and is bounded by the free surface of the water at $x_{3} = 0$ and the wetted surface of the twin hull. In the body-fixed coordinate system $(x_{1}, x_{2}, x_{3})$ , the forward motion of the vessel at speed $U$ along the $x_{1}$ -direction translates to the existence of a uniform water flow moving at constant speed $U$ opposite to this direction.

The dynamic behavior of the vessel is analyzed in the framework of linear wave theory, enabling the flow variables to be decomposed into a steady background flow and a time-dependent wave component, which is treated as harmonic. The velocity fields of the total flow components are described by the gradients of appropriate potential functions defined in $Ω$ .

On the basis of linear theory, each of these motions induces a distinct outgoing (radiated) wave field that contributes to the unsteady hydrodynamic responses. In this work’s context, $Φ$ is employed to denote potential functions of time-dependent flows, while (lowercase) $φ$ is utilized to represent the potential of time-independent flows, ensuring clear differentiation between the two. Furthermore, potential functions involved in the steady problem are annotated with the superscript $(S)$ , whereas the superscript $(U)$ is used for potential functions corresponding to the wave problem. This serves to amplify the distinction between the two problems, since time-independent (complex) potential functions are also used for the wave problem, exploiting the assumption of harmonic time dependence, as explained in greater detail in the sequel.

2.1. Formulation of the Steady Problem

The steady flow is defined by the interaction of a uniform parallel flow

U_{\infty} = (- U_{\infty}, 0, 0)

directed towards the negative

x_{1}

-axis, with the wetted surface in the domain

Ω

(see Figure 3). The perturbation field is calculated following a Neumann–Kelvin (NK) formulation (see, e.g., ref. [10]) using the decomposition,

$φ^{(S)} (x) = - U_{\infty} x_{1} + φ_{d}^{(S)} (x) .$

(1)

The potential function

φ_{d}^{(S)} (x)

is used to describe the generated Kelvin pattern and is expected to exhibit wave-like behavior downstream. The latter function is obtained as solution to the Laplace equation, supplemented by the free-surface boundary condition and an impermeability condition at the wetted surface. Considering the transom stern geometry of the studied vessel, a false body model is adopted (see e.g., ref. [11]) to ensure that the flow separates smoothly from the hull’s surface at the aft end. This approach is based on the addition of a virtual appendage (false body—FB) to the transom. The latter encloses the separated flow at low speeds, or the created air cavity in the higher speed range [12], and excludes this region from the water body simulated by the potential flow solver. The false body naturally extends the hull geometry, initiating tangentially from the entire stern profile, ensuring a smooth geometry and thereby maintaining the integrity of the flow at the aft end of the vessel. This approach presumes the transom remains fully ventilated at all speeds, leading to transom stern drag, attributed to the absence of hydrostatic pressure on the stern, to be overestimated at low speeds [13]. The false body length is set to

λ^{(S)} / 2

, where

λ^{(S)} = 2 π U_{\infty}^{2} / g

is the transverse wavelength of the Kelvin pattern. However, in cases where

λ^{(S)} / 2

exceeds

3 B_{S}

B_{S}

being the stern beam at mean water level (MWL), the virtual appendage length is limited to

3 B_{S}

. This restriction is based on experimental data provided in ref. [14], where it is shown that the boundary layer typically reattaches at a distance equivalent to six times the step height

(B_{S} / 2)

for high Reynolds number turbulent flows. Nonetheless, relevant calculations indicated that the evaluated wave-making resistance is not severely affected by this length, regardless of the precise value within these limits.

Considering the addition of the virtual boundary

\partial Ω_{F B}

, the potential function

φ_{d}^{(S)} (x)

is evaluated as a solution to the following BVP,

$\nabla^{2} φ_{d}^{(S)} (x) = 0, x \in Ω,$

(2)

$\frac{\partial^{2} φ_{d}^{(S)} (x)}{\partial {x_{1}}^{2}} + \frac{g}{U_{\infty}^{2}} \frac{\partial φ_{d}^{(S)} (x)}{\partial x_{3}} = 0, x \in \partial Ω_{F S},$

(3)

$\frac{\partial φ_{d}^{(S)} (x)}{\partial n} = - U_{\infty} n = U_{\infty} n_{1}, x \in \partial Ω_{W S} \cup \partial Ω_{F B},$

(4)

supplemented by appropriate conditions at infinity. In the above equations, $n = (n_{1}, n_{2}, n_{3})$ represents the unit vector, normal to the boundary $\partial Ω$ , pointing outward from the water body, as also shown in Figure 3. For the treatment of the above BVP, a low-order panel method is adopted using a simple singularity (source) distribution on $\partial Ω$ (see, e.g., ref. [15]), in conjunction with an integral representation of $φ_{d}^{(S)} (x)$ . Exploiting the inherent symmetry of the discussed (steady) problem with respect to the plane $x_{2} = 0$ , the following integral representation is utilized, which makes use of a mirroring technique to amplify computational efficiency,

$φ_{D}^{(S)} (x) = \int_{\partial \hat{Ω}} μ^{(S)} (x^{'}) G (x^{'}, x) d S (x^{'}), x \in \tilde{Ω}, x^{'} \in \partial \tilde{Ω},$

(5)

$where \tilde{Ω} = \{x \in Ω| x_{2} \geq 0\} and$

(6)

$G (x^{'}, x) = \frac{1}{4 π} (\frac{1}{|x^{'} - x|} + \frac{1}{|x^{'} - \tilde{x}|}) .$

(7)

The Green’s function of the Laplace equation in 3D, described by Equation (7), involves the mirror point of $x$ with respect to the symmetry plane: $\tilde{x} = (x_{1}, - x_{2}, x_{3}),$ and $μ^{(S)}$ is a source/sink strength distribution, defined on $\partial \tilde{Ω}$ . The above technique introduces an artificial boundary at the symmetry plane, where a homogeneous Neumann BC is intrinsically satisfied due to the Green’s function being used. The above formulation is combined with an appropriate scheme to satisfy the conditions at infinity based on the discrete Dawson operator, as described in the sequel.

2.2. Three-Dimensional BEM for the Steady-Flow Problem

The geometry of the different sections of

\partial \tilde{Ω}

is approximated using 4-node quadrilateral elements, on which the strength of the singularity distribution is considered piecewise constant. In the discrete model, the field equation is inherently satisfied by the superposition of the fundamental fields generated by all elements, while the boundary conditions are satisfied at the centroid of each panel (collocation points). The induced potential

φ_{k, j}

and velocities

(U_{k, j}, V_{k, j}, W_{k, j})

associated with the

j

-element’s contribution to the

k

-collocation point are numerically calculated, and the corresponding matrices of induced potential (

φ)

and velocity

(U, V, W)

, respectively, are computed. The latter square matrices have dimension

M = M_{F S} + M_{W S} + M_{F B}

, where

M_{F S}

M_{W S}

, and

M_{F B},

respectively, represent the number of quadrilateral elements distributed on the parts of

\partial \tilde{Ω},

corresponding to the free surface and the actual wetted surface and the false wetted surface, ensuring global continuity of geometry approximation. Consequently, the BVP is reduced to the algebraic system,

$\sum_{j = 1}^{M} A_{k, j}^{(S)} μ_{j}^{(S)} = b_{k}^{(S)}, k = 1, 2, \dots, M,$

(8)

whose solution provides the strength of the source/sink distribution $μ^{(S)}$ on each element. The latter, in conjunction with Equation (5), fully describes the resulting function $φ_{d}^{(S)} (x)$ in $\tilde{Ω}$ , while for $x_{2} < 0,$ it holds that $φ_{d}^{(S)} (x_{1}, x_{2}, x_{3}) = φ_{d}^{(S)} (x_{1}, - x_{2}, x_{3})$ . The numerical scheme involves a streamline-like arrangement of the boundary elements on $\partial {\tilde{Ω}}_{F S}$ , as shown in Figure 4. The acceleration in the $x_{1}$ -direction, which is involved in the free surface boundary condition [Equation (3)], is approximated by the derivative of the velocity in the $ξ$ -direction, as depicted in the details shown in Figure 4b,c.

$\frac{\partial^{2} φ_{d}^{(S)} (x)}{\partial x_{1}^{2}} = \frac{\partial U_{d}^{(S)} (x)}{\partial x_{1}} \approx - \frac{\partial U_{d}^{(S)} (x)}{\partial ξ} .$

(9)

The present model involves a four-point upstream finite difference (FD) scheme based on the Dawson backward operator [16], which is defined as follows:

${\frac{\partial U}{\partial ξ}|}_{k} \approx A_{k} \cdot U_{k} + B_{k} \cdot U_{k - 1} + C_{k} \cdot U_{k - 2} + D_{k} \cdot U_{k - 3},$

(10)

where

$\begin{array}{l} E_{k} = & (ξ_{k - 1} - ξ_{k}) \cdot (ξ_{k - 2} - ξ_{k}) \cdot (ξ_{k - 3} - ξ_{k}) \cdot (ξ_{k - 3} - ξ_{k - 1}) \cdot (ξ_{k - 2} - ξ_{k - 1}) \cdot \\ (ξ_{k - 3} - ξ_{k - 2}) \cdot (ξ_{k - 3} + ξ_{k - 2} + ξ_{k - 1} - 3 ξ_{k}), \end{array}$

(11)

$D_{k} = \frac{1}{E_{k}} [{(ξ_{k - 1} - ξ_{k})}^{2} \cdot {(ξ_{k - 2} - ξ_{k})}^{2} \cdot (ξ_{k - 2} - ξ_{k - 1}) \cdot (ξ_{k - 2} + ξ_{k - 1} - 2 ξ_{k})],$

(12)

$C_{k} = - \frac{1}{E_{k}} [{(ξ_{k - 1} - ξ_{k})}^{2} \cdot {(ξ_{k - 3} - ξ_{k})}^{2} \cdot (ξ_{k - 3} - ξ_{k - 1}) \cdot (ξ_{k - 3} + ξ_{k - 1} - 2 ξ_{k})],$

(13)

$B_{k} = \frac{1}{E_{k}} [{(ξ_{k - 2} - ξ_{k})}^{2} \cdot {(ξ_{k - 3} - ξ_{k})}^{2} \cdot (ξ_{k - 3} - ξ_{k - 2}) \cdot (ξ_{k - 3} + ξ_{k - 2} - 2 ξ_{k})],$

(14)

$A_{k} = - (B_{k} + C_{k} + D_{k}) .$

(15)

In the above equations,

ξ

denotes the distance covered in a curve defined in the

ξ

-direction, measured from an arbitrary point, and the velocities involved in Equation (10) are equal to

$U_{k} = \sum_{j = 1}^{M} μ_{j}^{(S)} \cdot U_{k, j} .$

(16)

Figure 4 depicts an indicative sparse boundary mesh, where it can be seen that additional grid lines are introduced behind the stern to minimize the deviation between the differentiation direction and the actual

x_{1}

-direction. Semicosine spacing is employed to discretize the wetted surface in the

x_{3}

-direction to achieve increased grid resolution at the keel, where the geometry exhibits the highest gradient. Furthermore, ghost nodes are employed to assess the induced acceleration at the first three element rows of the grid (see Figure 4b). The value of the induced velocity at the ghost nodes is taken to be equal to the velocity induced at the centroid of the last actual element involved in the differentiation. The same technique is also applied to compute the acceleration at the first three element rows behind the false body mesh.

Based on the above finite difference scheme, the (non-square) matrix

\partial U

of dimensions

M_{F S} \times M

is defined. The elements of this matrix,

$\partial U_{k, j} = A_{k} \cdot U_{k, j} + B_{k} \cdot U_{k - 1, j} + C_{k} \cdot U_{k - 2, j} + D_{k} \cdot U_{k - 3, j}$

(17)

model the contribution of the $j$ -element on $\partial U_{k} / \partial ξ$ at the $k$ -collocation point on the free surface. Consequently, the discrete model of the steady-flow problem, modeling the disturbance field, is defined by the linear algebraic system of Equation (8), where

$A_{k, j}^{(S)} = \{\begin{matrix} \partial U_{k, j} + \frac{g}{U_{\infty}^{2}} W_{k, j}, & E l e m e n t (k) \in \partial Ω_{F S} \\ U_{k, j} n_{k, 1} + V_{k, j} n_{k, 2} + W_{k, j} n_{k, 3}, & E l e m e n t (k) \in (\partial Ω_{W S} \cup \partial Ω_{F B}) \end{matrix},$

(18)

$b_{k}^{(S)} = \begin{matrix} 0, & E l e m e n t (k) \in \partial Ω_{F S} \\ U_{\infty} n_{k, 1}, & E l e m e n t (k) \in \partial (\partial Ω_{W S} \cup \partial Ω_{F B}) . \end{matrix}$

(19)

After obtaining the solution to the above linear system, the free surface elevation can be obtained from the dynamic boundary condition as

$η_{d}^{(S)} (x) = \frac{U_{\infty}}{g} \frac{\partial φ_{d}^{(S)} (x_{1}, x_{2}, 0)}{\partial x_{1}}, x \in \partial D_{F S} .$

(20)

An indicative result of the field generated by the vessel moving at 10 knots is shown in Figure 5. Specifically, Figure 5a illustrates the wave-like behavior of the steady perturbation field trailing behind the vessel, demonstrating the characteristic wake pattern, and Figure 5b provides a detailed view of the solution on the computational grid in a localized area around the ship, also highlighting the directly computed solution and its mirrored counterpart. The resulting dynamic pressure on the wetted surface for the same velocity is shown in Figure 6, illustrating the pressure distribution along the hull.

The results refer to a draft of

1.55 m

, ensuring the availability of experimental data for verification purposes. The latter data concerning calm water resistance can be found in ref. [9].

The total calm water resistance is approximated by the sum of the wave-making resistance and the frictional resistance,

$R_{T} = R_{W} + R_{F} = [C_{W} + (1 + β k) C_{F}] \frac{1}{2} ρ U_{\infty}^{2} A_{W S}$

(21)

where $C_{W}$ and $C_{F},$ respectively, denote the wave-making and frictional resistance coefficients, $(1 + β k)$ is a form factor and $A_{W S}$ is the wetted surface area. The coefficient $C_{W}$ is computed using data by the steady BEM model discussed earlier. In particular, $C_{W}$ equals

$C_{W} = \frac{1}{A_{W S}} \int_{\partial Ω_{W S}} \frac{p - p_{\infty}}{0.5 ρ U_{\infty}^{2}} \cdot n_{1} d S = \frac{1}{A_{W S}} \int_{\partial Ω_{W S}} (1 - \frac{V^{2} (x) - 2 g z}{U_{\infty}^{2}}) \cdot n_{1} d S$

(22)

where

$V (x) = {[{(- U_{\infty} + \frac{\partial φ_{d}^{(S)} (x)}{\partial x_{1}})}^{2} + {(\frac{\partial φ_{d}^{(S)} (x)}{\partial x_{2}})}^{2} + {(\frac{\partial φ_{d}^{(S)} (x)}{\partial x_{2}})}^{2}]}^{1 / 2} .$

(23)

The frictional resistance coefficient

C_{F}

is evaluated using the ITTC 1957 formula [8],

$C_{F} = \frac{0.075}{{(\log (Re) - 2)}^{2}}$

(24)

and is scaled by the form factor $(1 + β k),$ which depends on the hull form and, in the case of a twin-hull configuration, can be approximated by the following equation based on the demi-hull’s slenderness ratio [17],

$1 + β k = 3.03 {(\frac{L}{\nabla^{1 / 3}})}^{- 0.4},$

(25)

where $\nabla$ denotes the submerged volume of a single demi-hull. For the considered twin-hull geometry and draft the form factor equals $(1 + β k) = 1.39$ . The Reynolds number is evaluated, considering the dynamic viscosity of seawater at temperature T = 20 °C, $(μ = 1.09 \times 10^{- 3} {Nsm}^{- 2})$ .

Figure 7 depicts the calm water resistance components along with the total calm water resistance, as evaluated by the present BEM scheme, and the results are compared against experimental data. As it can be observed in Figure 7, the wave-making resistance is augmented by a stern drag term resulting from the absence of hydrostatic pressure at the aft end, leading to an overestimation of the total resistance at low Froude numbers, as expected from the present scheme. In the context of linearized theory, this drag remains unchanged at different speeds. Apart from that, the results demonstrate good agreement with the experimental data within the range of Froude number

F n < 0.25

. For higher

F n

values, the present model underestimates the resistance. However, significantly better agreement with the experimental data can be obtained by considering the nonlinear problem, including dynamic trim effects. These findings will be further explored in future extensions of the present work.

2.3. Formulation of the Unsteady Problem

The dynamic motions of the vessel, based on the prevailing sea conditions, are derived using standard linear hydrodynamic analysis in the frequency domain. In particular, assuming that all time-dependent quantities oscillate harmonically in the form

\exp (i ω_{e} t

), where

ω_{e}

denotes the encounter frequency in

[rad / s]

and

i = \sqrt{- 1}

, the problem is transferred to the frequency domain using the following representation for the total unsteady potential function,

$\begin{matrix} Φ^{(U)} (x; β, t) = Φ_{0}^{(U)} (x; β, U_{\infty}, t) + Φ_{d}^{(U)} (x; β, U_{\infty}, t) + \sum_{m = 1}^{6} \frac{d ξ_{m}}{d t} Φ_{m}^{(U)} (x; β, U_{\infty}, t) = \\ Re \frac{i g H}{2 ω_{0}} φ^{(U)} (x; β, U_{\infty}) \cdot \exp (i ω_{e} t) . \end{matrix}$

(26)

In Equation (26),

H = 2 A

denotes the incident wave height,

g

is the gravitational acceleration,

ξ_{m}

is the complex amplitude of the vessel’s response towards the

m

-th generalized direction,

β

is the angle of incidence of the wave field with respect to the

x_{1}

-axis and

U_{\infty}

is the vessel’s velocity in the

x_{1}

-direction. The encounter frequency

ω_{e}

is defined as

$ω_{e} = \{{φ_{0}}^{(U)} (x; β, U_{\infty}) + {φ_{d}}^{(U)} (x; β, U_{\infty}) + \sum_{m = 1}^{6} ξ_{m} {φ_{m}}^{(U)} (x; β, U_{\infty}) .$

(28)

In Equation (28), $i ω_{e} φ_{0}^{(U)} (x, β, U_{\infty})$ and $i ω_{e} φ_{d}^{(U)} (x, β, U_{\infty})$ , respectively, stand for complex amplitudes of the incident and the diffracted subfields. Furthermore, $φ_{m}^{(U)} (x, β, U_{\infty}), m = 1, 2, \dots, 6$ is the complex potential of the radiation field generated by a unit amplitude oscillation of the vessel along the $m^{t h}$ generalized direction.

The incident potential function

φ_{0}^{(U)} (x, β, U_{\infty})

is considered known and equal to

${φ_{0}}^{(U)} (x; β, U_{\infty}) = \frac{g A}{ω_{0} ω_{e}} \exp (k_{0} x_{3}) \exp [- i k_{0} (\cos (β) x_{1} + \sin (β) x_{2}] .$

(29)

The unsteady field is computed by modeling the whole domain, since it does not present symmetries, except in the special cases where

β

equals

0

180

degrees. Furthermore, the transom sterns of both demi-hulls are modeled as vertical boundaries, incorporated in

\partial Ω_{W S}

, omitting the false body representation discussed in the previous section. The diffraction and the six radiated fields are obtained as solutions to the following BVPs:

$\nabla^{2} {φ_{k}}^{(U)} (x) = 0, k = d, 1, 2, \dots, 6, x \in Ω,$

(30)

$\begin{array}{l} \frac{\partial {φ_{k}}^{(U)} (x)}{\partial x_{1}} - μ (x) {φ_{k}}^{(U)} (x) + \frac{{U_{\infty}}^{2}}{g} \frac{\partial^{2} {φ_{k}}^{(U)} (x)}{\partial {x_{1}}^{2}} \\ - \frac{2 i ω_{e} U_{\infty}}{g} \frac{\partial {φ_{k}}^{(U)} (x)}{\partial x_{1}} = 0, k = d, 1, 2, \dots, 6, x \in \partial Ω_{F S}, \end{array}$

(31)

$\frac{\partial {φ_{d}}^{(U)} (x)}{\partial n} = \frac{\partial {φ_{0}}^{(U)} (x)}{\partial n}, x \in \partial Ω_{W S},$

(32)

$\frac{\partial {φ_{k}}^{(U)} (x)}{\partial n} = N_{k} + M_{k}, k = 1, 2, \dots, 6, x \in \partial Ω_{W S} .$

(33)

In Equation (33),

N_{k}

denotes the components of the generalized normal vector, defined as

N_{k} = n_{k}, N_{k + 3} = {(n \times x)}_{k}, k = 1, 2, 3

. Consequently, the vessel’s

k^{t h}

DoF is defined as translational motion along the

k^{t h}

axis for

k = 1, 2, 3

or rotational motion around the axis

(k - 3)

for

k = 4, 5, 6

. The quantities

M_{k}

, appearing in Equation (33) are defined as functions of the derivatives of the relative flow velocity of the steady field at the mean position of the wetted surface. The latter quantities mainly involve higher-order derivatives of the perturbation potential due to the motion of the ship at a constant speed (see e.g., ref. [18]), and after linearization, they can be approximated by

$M_{k} \approx 0, k = 1, 2, 3, 4, M_{5} \approx \frac{U_{\infty}}{i ω_{e}} n_{3}, M_{6} \approx - \frac{U_{\infty}}{i ω_{e}} n_{2} .$

(34)

Numerical solutions to the above BVPs are obtained using a low-order BEM, which employs piecewise constant source distributions on quadrilateral boundary elements. The seven unknown fields are represented by the following integral:

$φ_{k}^{(U)} (x) = \int_{\partial \hat{Ω}} {μ_{k}}^{(U)} (x^{'}) G (x^{'}, x) d S (x^{'}), x \in Ω, x^{'} \in \partial Ω, k = d, 1, 2, \dots, 6, where$

(35)

$G (x^{'}, x) = {|x^{'} - x|}^{- 1} / 4 π .$

(36)

The above BVPs are treated by a low-order panel method based on simple singularity (source) distribution on quadrilateral elements. The radiating behavior of the diffraction, as well as the six radiation subfields at infinity, is treated by adopting a Perfectly Matched Layer (PML) technique. This technique is utilized to attenuate the evaluated solutions far from the vessel’s position and prevent the occurrence of numerical reflections from the outer layers of the computational mesh. The PML is implemented by complexifying the frequency parameter

μ = ω_{e}^{2} / g

, involved in the Free Surface BC [Equation (31)], using an appropriately defined imaginary component. The parameter is redefined as

$μ (x; ω_{e}) = \{\begin{matrix} {ω_{e}}^{2} g^{- 1}, |x| < R_{P M L} (θ (x_{1}, x_{2})) \\ {ω_{e}}^{2} g^{- 1} {(1 + i c \frac{{[|x| - R_{P M L} (θ (x_{1}, x_{2}))]}^{n}}{λ_{e}^{n}})}^{2}, |x| \geq R_{P M L} (θ (x_{1}, x_{2})) . \end{matrix}$

(37)

In Equation (37),

λ_{e}

stands for the wavelength that propagates at frequency

ω_{e}

and is not affected by the background flow. This corresponds to the wavelength of the diffraction and radiation fields along the

x_{2}

-direction. Conversely, the wavelengths of the field components propagating along the

x_{1}

-direction are significantly affected by the background flow. Specifically, the wavelengths aft of the stern are elongated, while compressed wavelengths are generated ahead of the bow in the subcritical case (

ω_{0} U_{\infty} / g < 1 / 4

for head seas). In supercritical cases, the diffraction and radiation subfields do not include components propagating toward the positive

x_{1}

-axis. This can be clearly observed in Figure 8, which illustrates the real part of the free surface elevation generated by the unit amplitude heaving oscillation of the vessel for a fixed absolute frequency corresponding to

λ_{0} / L = 1

, at two different speeds, respectively falling into the subcritical (a) and supercritical (b) regimes. Optimal values for the parameters c and n, involved in Equation (37), are selected for minimizing the numerical reflections; see [19]. The PML activation curves are also depicted in Figure 8 using dashed lines. The evaluation of the acceleration in the

x_{1}

-direction, which is involved in Equation (31), is achieved by using the discrete Dawson operator, described in the previous section. This suggests that the boundary mesh used for evaluating the involved unsteady subfields also comprises streamline-like lines on the free surface resulting in a rectangular mesh geometry.

After obtaining all the subfields, the response vector

(ξ)

of the vessel in 6 DoF is evaluated as the solution to the following system of equations:

$[- {ω_{e}}^{2} (M + A (ω_{e})) + i ω_{e} (B (ω_{e}) + U_{\infty} N) + C] ξ = F_{0} + F_{d} .$

(38)

In Equation (38), the inherent inertia of the vessel is modeled via the

6 \times 6

tensor

M

, which equals

$M = (\begin{matrix} M^{t t} & M^{t r} \\ M^{r t} & M^{r r} \end{matrix}),$

(39)

where the diagonal matrix $M^{t t}$ models the inertia of the translational degrees of freedom with respect to the applied forces and is defined as

$M_{i i}^{t t} = M = ρ V = ρ \int_{\partial Ω_{W S}} x_{3} n d S, i = 1, 2, 3,$

(40)

with $M$ denoting the vessel’s mass and $V$ the submerged volume (of both demihulls). The diagonal elements of the matrix $M^{r r}$ , which models the inertia of the rotational degrees of freedom with respect to the applied moments, equal the corresponding moments of inertia of the vessel with respect to the selected coordinate system and are evaluated as

$M_{i i}^{r r} = I_{i i} = M R_{i}^{2}, i = 1, 2, 3,$

(41)

where $R_{i}$ represents the radius of gyration about the $i^{t h}$ axis. In the following, the radii of gyration are defined as $R_{1} = 0.2 B$ , $R_{i} = 0.25 L, i = 2, 3,$ and the products of inertia $I_{12}$ and $I_{13}$ are assumed to be zero due to symmetry relative to the longitudinal axis. Additionally, due to lack of detailed information and assuming symmetry of the mass distribution relative to the beam axis $(x_{2})$ —even without explicit symmetry in the geometry—the product $I_{23}$ is also considered negligible, making the matrix $M^{r r}$ also diagonal. The submatrices $M^{t r}$ and $M^{t r}$ are defined as follows:

$M^{t r} = - M^{r t} = (\begin{matrix} 0 & J_{3} & - J_{2} \\ - J_{3} & 0 & J_{1} \\ J_{2} & - J_{1} & 0 \end{matrix}),$

(42)

where $J = \vec{O G} = (J_{1}, J_{2}, J_{3})$ is the position vector of the center of gravity $G$ with respect to the origin, which is taken to be the center of flotation. The longitudinal position of $G$ is taken to coincide with that of the center of buoyancy, which is evaluated as the center of the submerged volume, to ensure that the vessel maintains an even keel. The vertical position of the center of gravity is placed at one-third of the hull depth measured from the keel, since this positioning represents a common approximation in ship design that supports accurate stability.

The added inertia of the vessel is represented by the

6 \times 6

matrix

A

, where the elements

A_{k, l} (ω_{e})

correspond to the components of the radiation loads, induced by the

l^{t h}

field on the

k^{t h}

DoF, in phase with the acceleration. The elements

B_{k, l} (ω_{e})

of the

6 \times 6

hydrodynamic damping matrix

B

are determined by the components that are in phase with the velocity. The above radiation forces and moments are computed by integrating the pressure exerted by the

l^{t h}

radiation field on

\partial Ω_{W S}

, multiplied by the component of the normal vector corresponding to the

k^{t h}

DoF. Therefore, the hydrodynamic coefficients (added mass and hydrodynamic damping) are computed using the following equation

$A_{k l} + \frac{B_{k l}}{i ω_{e}} = ρ \int_{\partial Ω_{W S}} (φ_{l}^{(U)} (x) - \frac{U_{\infty}}{i ω_{e}} \frac{\partial φ_{l}^{(U)} (x)}{\partial x_{1}}) n_{k} d S (x), k, l = 1, 2, \dots, 6 .$

(43)

It is noted that the above matrices are dependent on the encounter frequency, and therefore, they also depend on the direction of propagation of the incident field, apart from the absolute frequency and the speed at which the vessel moves. The hydrostatic restoring forces are modeled via Matrix

C

. As concerns the transitional DoFs, only vertical movements induce hydrostatic restoration and thus

C_{33} > 0

. Regarding the rotational DoFs, since the origin has been selected to coincide with the center of floatation (centroid of the water plane), all the first-order moments as well as the products of inertia of the waterplane are zero. Based on the above observations, the nonzero elements of

C

are the following:

$C_{33} = ρ g A_{W P},$

(44)

$C_{44} = M g \bar{G B} + ρ g \int_{W P} {x_{2}}^{2} d S = M g \bar{G M_{1}},$

(45)

$C_{55} = M g \bar{G B} + ρ g \int_{W P} {x_{1}}^{2} d S = M g \bar{G M_{2}},$

(46)

where $W P$ denotes the waterplane, $A_{W P}$ is the waterplane surface, $\bar{G B}$ is the vertical distance between the center of buoyance and the center of gravity and $G M_{k}, k = 1, 2$ is the metacentric height of the vessel with respect to the $k^{t h}$ axis. The additional $6 \times 6$ matrix $N$ , defined as,

$N = (\begin{matrix} 0 & N^{t r} \\ 0 & N^{r r} \end{matrix}), where N^{t r} = (\begin{matrix} 0 & 0 & 0 \\ 0 & 0 & M \\ 0 & - M & 0 \end{matrix}), N^{r r} = (\begin{matrix} 0 & - J_{2} & - J_{3} \\ 0 & J_{1} & 0 \\ 0 & 0 & J_{1} \end{matrix}),$

(47)

represents additional damping and coupling effects on the dynamic behavior of the vessel due to forward motion. In the context of linear theory, the effect of $N$ arises since products of the form $U_{\infty} ξ_{k}$ are comparable to linear terms in the equations of motion and thus cannot be neglected after linearization; see [18,20,21].

Finally, the forces and moments acting on the hull are evaluated by integration of the pressure induced by the incident and the diffracted subfields on the wetted surface. In particular, the components of the Froude–Krylov and diffraction forces are computed, taking into account the vessel’s forward motion by the following equation [20,21]:

$F_{k l} = ρ ω_{e}^{2} \int_{\partial Ω_{W S}} (φ_{l}^{(U)} (x) - \frac{U_{\infty}}{i ω_{e}} \frac{\partial φ_{l}^{(U)} (x)}{\partial x_{1}}) n_{k} d S (x), k = 1, 2, \dots, 6, l = 0, d .$

(48)

Source link

Alexandros Magkouris www.mdpi.com

Greenberg News

A Novel BEM for the Hydrodynamic Analysis of Twin-Hull Vessels with Application to Solar Ships

2.1. Formulation of the Steady Problem

2.2. Three-Dimensional BEM for the Steady-Flow Problem

2.3. Formulation of the Unsteady Problem

Greenberg

2.1. Formulation of the Steady Problem

2.2. Three-Dimensional BEM for the Steady-Flow Problem

2.3. Formulation of the Unsteady Problem

Related Posts

Melatonin Improves H2O2-Induced Oxidative Stress in Sertoli Cells Through Nrf2-Keap1 Signaling Pathway

PEPH Newsletter | National Institute of Environmental Health Sciences

Environmental Factor – May 2023: Reducing phthalates in beauty products may lower health risks, disparities

Greenberg