Abstract

Translate

A variety of Offshore Floating Photovoltaics (OFPVs) applications rely on the capacity of their floating support structures displacing in the shape of surface waves to reduce extreme wave-induced loads exerted on their floating-mooring system. This wave-adaptive displacement behaviour is typically realized through two principal design approaches, either by employing slender and continuously deformable structures composed of highly elastic materials or by decomposing the structure into multiple floating rigid pontoons interconnected via flexible connectors. The hydrodynamic behaviour of these structures is commonly analyzed in the literature using potential flow theory, to characterize wave loading, whereas in order to deploy such OFPV prototypes in realistic marine environments, a high-fidelity numerical fluid–structure interaction model is required. Thus, a versatile three-dimensional numerical scheme is herein presented that is capable of handling non-linear fluid-flexible structure interactions for Very Flexible Floating Structures (VFFSs): Multibody Dynamics (MBD) for modularized floating structures and floating-mooring line interactions. In the present study, this is achieved by employing the Smoothed Particles Hydrodynamics (SPH) fluid model of DualSPHysics, coupled both with the MBD module of Project Chrono and the MoorDyn+ lumped-mass mooring model. The SPH-MBD coupling enables modelling of large and geometrically non-linear displacements of VFFS within an Applied Element Method (AEM) plate formulation, as well as rigid body dynamics of modularized configurations. Meanwhile, the SPH-MoorDyn+ captures the fully coupled three-dimensional response of floating-mooring and floating-floating dynamics, as it is employed to model both moorings and flexible interconnectors between bodies. The coupled SPH-based numerical scheme is herein validated against physical experiments, capturing the hydroelastic response of VFFS, rigid body hydrodynamics, mooring line dynamics, and flexible connector behaviour under wave loading. The demonstrated numerical methodology represents the first validated Computational Fluid Dynamics (CFD) application of moored VFFS in three-dimensional domains, while its robustness is further confirmed using modular floating systems, enabling OFPV engineers to comparatively assess these two types of wave-adaptive designs in a unified numerical framework.

Full text

Turn on search term navigation

Translate

1. Introduction

Efforts to commercially deploy offshore solar technologies are still ongoing, as standardized design frameworks have yet to be established. Unlike the more mature offshore wind technologies, current prototypes of Offshore Floating Photovoltaics (OFPVs) exhibit significant diversity, particularly in the types of their floating support structures [1,2,3]. One prominent design criterion in some of the developed frameworks focuses on leveraging the displacement behaviour of the floating support structure, enabling it to almost fully follow the motion of surface waves. This dynamic response promotes a more uniform force distribution along the structure’s length and significantly reduces the maximum wave-induced forces experienced by its mooring system.

A widely adopted approach to introduce such displacement behaviour in floating support structures involves the use of rigid pontoons interconnected via flexible connectors, collectively forming a semi-rigid deformable system that supports solar panels [4,5,6]. The individual pontoons, lightweight and in direct contact with the water surface, are interconnected using linear soft [7,8], linear stiff [9], or hinged connectors [10,11,12], depending on the desired mechanical response. The hydrodynamic performance of such modular arrays has been extensively studied within the framework of linear potential flow theory [13,14,15,16,17,18], which limits their applicability under highly non-linear offshore wave conditions. Time-domain analyses have also been conducted using linear [19] and higher-order Stokes wave theories [8], yet these still preclude the accurate representation of non-linear Fluid–Structure Interaction (FSI) phenomena. In contrast, Computational Fluid Dynamics (CFD) models have been successfully applied to capture complex non-linear fluid–multibody interactions in elevated modular OFPV configurations [20,21], as well as in wave energy converter arrays [22], shared mooring systems involving multiple floating bodies [23], and moored single-body offshore wind platforms [24,25]. Nevertheless, application of CFD models in modularized OFPV systems resting on the water surface is scarce, with many studies relying on potential flow models enhanced by CFD tuned [26,27] or experimental tuned [16,18] viscous damping parameters, to allow for computationally efficient numerical investigations. Thus, the present study aims to develop a computationally feasible methodology for fully CFD resolved wave-deformable modularized OFPV structures by employing a variable resolution CFD scheme for the descritized fluid domain [28,29].

Alternatively, very thin floating plates, in combination with very elastic materials, utilize their inherent structural flexibility to deform relative to wave surface elevations for OFPV applications [30,31]. These systems, referred to as Very Flexible Floating Structures (VFFSs) [32], can be analyzed using Computational Fluid Dynamics (CFD) models formulated in either Eulerian [33,34,35,36,37] or Lagrangian frameworks [38,39,40,41]. The latter, due to their inherent ability to deal with floating boundaries subjected to large displacements, are particularly suited for these systems. Among particle-based schemes, Smoothed Particle Hydrodynamics (SPH) has become the most widely adopted approach [42], enabling the representation of structural elasticity through both monolithic [43,44,45,46] and partitioned [40,41,47,48,49,50,51] schemes. However, applications of SPH to thin flexible floating structures have largely remained two-dimensional, with fully three-dimensional applications demonstrated only recently in SPH [41]. However, VFFS in fully three-dimensional moored configurations are still lacking for both Lagrangian and Eulerian methods, underscoring a research gap for CFD applications, which mostly focus on moored rigid floating structures [52].

To enable the modelling of all the above applications, a versatile three-dimensional numerical framework is developed, capable of resolving the following components:

(a). Non-linear fluid dynamics, addressed through the SPH solver of DualSPHysics [53], where a variable resolution scheme is employed to efficiently extend the model to fully three-dimensional domains while maintaining numerical accuracy [28,29].

(b). Multibody Dynamics (MBD), handled through the multiphysics library Project Chrono [54], which solves the rigid-body dynamics of modular components. This module is coupled with the SPH solver to capture non-linear fluid–structure interactions of multibody systems [55].

(c). Structural flexibility of thin floating plates, represented using the Applied Element Method (AEM) incorporated within the SPH–MBD coupling [40,41], enabling the simulation of large structural displacements. The implemented AEM formulation has been validated for both beam and plate configurations under static, dynamic, and hydrodynamic conditions.

(d). Mooring line dynamics, resolved through the MoorDyn+ lumped-mass model [56], an enhanced implementation of the original MoorDyn framework [57]. The model resolves the dynamic equilibrium of mooring lines and flexible connectors under external hydrodynamic forcing while being two-way coupled with the SPH solver [52] to accurately capture interactions between moored floating structures and non-linear wave excitation.

The coupled SPH-based numerical framework is validated both quantitatively and qualitatively against physical experiments involving two representative configurations, an array of interconnected rigid pontoons, and a uniform flexible floating structure. Through these validations, the model demonstrates its capability to accurately capture the three-dimensional hydroelastic response under wave loading of flexible and modularized floating structures along with the corresponding mooring line tensions.

2. Numerical Methods

The equations of the SPH fluid model employed in this study are presented in Section 2.1.1, while the description of the domain behaviour within the variable resolution scheme is provided in Section 2.1.2. The rigid body dynamics equations and the corresponding mechanical constraints are detailed in Section 2.2.1, along with their coupling procedure with the SPH solver (Section 2.2.2). The lumped-mass mooring model equations are presented in Section 2.3.1 and their coupling procedure with the SPH equations is described in Section 2.3.2. The formulation of the AEM plate within the SPH–MBD coupling is detailed in Section 2.4.

2.1. SPH Fluid Solver

2.1.1. Governing Equations

SPH equations, implemented in DualSPHysics [53] and solved using a Graphics Processing Unit (GPU), are outlined as adopted in this study. The fluid domain is discretized into a discrete set of particles, each representing a point at which physical properties are assigned, namely position $r$ , velocity $v$ , and density $ρ$ . These particles are initially spatially distributed with a uniform inter-particle spacing, denoted as dp, throughout the computational domain. Their associated properties evolve over time according to the Lagrangian form of the Navier–Stokes equations, which govern the conservation of momentum (Equation (1)) and mass (Equation (2)). In these equations, each fluid particle a interacts with $N_{b}$ surrounding fluid particles b within a kernel function $W_{a b}$ [58], defined by a smoothing length $h = 1.5 d p$ and cut-off radius $2 h$ . Particle pressures P are evaluated according to particle densities using an equation of state, presented in Equation (3) [59], where the isotropic factor is assigned a value of 1.

(1) $\begin{matrix} \frac{d v_{a}}{d t} = \underset{Non - Conservative Pressure Term}{\underset{⏟}{- \sum_{b = 1}^{N_{b}} m_{b} (\frac{f (P_{a}, P_{b})}{ρ_{a} ρ_{b}}) \nabla_{a} W_{a b}}} + \underset{Gravity Term}{\underset{⏟}{g}} + \underset{Artificial Viscosity Term}{\underset{⏟}{Π_{a}}} + \dots \\ \dots + \underset{Arbitrary Lagrangian - Eulerian Formulation}{\underset{⏟}{\sum_{b = 1}^{N_{b}} \frac{m_{b}}{ρ_{b}} (v_{a} \otimes δ v_{a} + v_{b} \otimes δ v_{b} + v_{a} \otimes (δ v_{a} - δ v_{b})) \cdot \nabla_{a} W_{a b}}} \end{matrix}$

(2) $\begin{matrix} \frac{d ρ_{a}}{d t} = \underset{Arbitrary Lagrangian - Eulerian Formulation}{\underset{⏟}{\sum_{b = 1}^{N_{b}} [\frac{ρ_{a}}{ρ_{b}} m_{b} (v_{a b} + δ v_{a b}) \cdot \nabla_{a} W_{a b} + \frac{m_{b}}{ρ_{b}} (ρ_{a} δ v_{a} + ρ_{b} δ v_{b}) \cdot \nabla_{a} W_{a b}]}} \\ + \underset{Density Diffusion Term}{\underset{⏟}{2 δ_{D} h c_{s} \sum_{b = 1}^{N_{b}} (ρ_{b a}^{T} - ρ_{a b}^{H}) \frac{r_{a b} \cdot \nabla_{a} W_{a b}}{r_{a b}^{2}} \frac{m_{b}}{ρ_{b}}}} \end{matrix}$

(3) $P = c_{s}^{2} (ρ - ρ_{0})$

where m is the particle mass,

g

the gravitational acceleration vector,

c_{s} = \sqrt{\partial P / \partial ρ}

the speed of sound, and

ρ_{0} = 1000 {kg / m}^{3}

the reference fluid density, while superscripts T and H refer to the total and hydrostatic components. In the mass equation, a density diffusion term is utilized, where the density diffusion coefficient

δ_{D}

is set equal to a default value of

0.10

based on previously conducted analytical and experimental validation cases [60], smoothing the density field and attenuating high-frequency pressure fluctuations.

In the momentum equation, a non-conservative pressure term (Equation (4)) is applied [61,62]. Particles located near or on the free surface region ( $Ω_{F} \cup Ω_{N F}$ ) follow a standard conservative pressure gradient. In contrast, inner fluid particles ( $Ω_{I}$ ), which do not interact with any free surface neighbouring particles, follow a linearized pressure gradient corrected by a factor $L_{a}$ (Equation (5)), suppressing tensile instability [63]. Free surface particles are detected using the approach of [64], while their normal vectors (Equation (6)) are computed through the particle shifting technique of [65], based on the concentration field gradient $\nabla C$ evaluated according to Equation (7). The computed normals are further smoothed by averaging over neighboring free surface particles normals.

(4) $f (P_{a}, P_{b}) = \{\begin{matrix} P_{a} + P_{b}, & if a \in Ω_{F} \cup Ω_{N F}, \\ (P_{b} - P_{a}) L_{a}, & if a \in Ω_{I} \end{matrix}$

(5) $L_{a} = {(\sum_{b = 1}^{N_{b}} \frac{m_{b}}{ρ_{b}} r_{b a} \otimes \nabla_{a} W_{a b})}^{- 1}$

(6) ${\bar{n}}_{a} = - \frac{L_{a} \cdot \nabla C_{a}}{∥ L_{a} \cdot \nabla C_{a} ∥}$

(7) $\nabla C_{a} = \sum_{b} \frac{m_{b}}{ρ_{b}} \nabla_{a} W_{a b}$

where

\bar{{(\cdot)}_{a b}} = [{(\cdot)}_{a} + {(\cdot)}_{b}] / 2

is the arithmetic mean and

{(\cdot)}_{a b} = [{(\cdot)}_{a} - {(\cdot)}_{b}]

is the arithmetic difference between particles a and b. Additionally, an artificial dissipation term

Π_{a}

is incorporated to model fluid viscous effects (Equation (8)), where the artificial viscosity coefficient

α

is set equal to a default value of

0.01

based on previously conducted experimental validation [66], and

η_{h}^{2} = 0.01 h^{2}

is a negligible non-zero term to prevent numerical singularities.

(8) $Π_{a} = \{\begin{cases} \sum_{b = 1}^{N_{b}} m_{b} (\frac{α c_{s}}{\bar{ρ_{a b}}}) (\frac{h v_{a b} \cdot r_{a b}}{r_{a b}^{2} + η_{h}^{2}}) \nabla_{a} W_{a b}, & v_{a b} \cdot r_{a b} < 0 \\ 0, & v_{a b} \cdot r_{a b} \geq 0 \end{cases}$

In both momentum and mass equations, particle velocities follow an Arbitrary Lagrangian–Eulerian (ALE) formulation [67], where corrections are applied to the pure Lagrangian velocities based on particle shifting techniques to ensure volume conservation. Position ( $δ r_{a}$ ) and velocity ( $δ r_{a}$ ) corrections (Equation (9)) are determined based on the shifting vector proposed by [68]. Subsequently, position shifting of particles is also applied, as shown in Equation (10).

(9) $δ r_{a} = δ v_{a} Δ t_{S P H} = D_{s} \nabla C_{a}$

(10) $\frac{d r_{a}}{d t} = v_{a} + δ v_{a}$

where

Δ t_{S P H}

the numerical timestep and

D_{s} = - 0.01 h^{2}

is a shifting diffusion coefficient.

The model employs modified dynamic boundary conditions based on the enhanced formulation proposed in [69]. Wave generation follows a second-order Stokes regular approach, implemented through a piston-type moving boundary whose motion is prescribed by AwaSys [70]. To mitigate wave reflection and energy accumulation, damping regions are introduced following the methodology described in [71]. Temporal advancement of the solution is carried out with a symplectic position–Verlet integration scheme [72]. The timestep varies dynamically according to the formulation of [73], ensuring compliance with the Courant–Friedrichs–Lewy ( $C F L$ ) criterion, which in this study is maintained at $C F L = 0.2$ .

2.1.2. Variable Resolution

A variable resolution approach is employed in this study to efficiently handle fluid–fluid and fluid–solid interactions, thereby reducing the computational cost associated with scaling to three-dimensional domains. The method, originally developed by [28,29], divides the global computational domain $Γ$ into a finite set of N subdomains $Γ_{i}$ , expressed as $Γ = ⋃_{i = 1}^{N} Γ_{i}$ , where $i = 1, \dots, N$ . Each subdomain satisfies the same governing equations presented in Section 2.1.1 but is discretized using its own initial inter-particle spacing $d p_{i} = 2 d p_{i + 1}$ and corresponding smoothing length $h_{i} = 1.5 d p_{i}$ . Adjacent subdomains overlap through buffer regions of thickness $2 h_{i}$ , following the open-boundary formulation proposed in [74]. These regions are populated with buffer particles that acquire their physical properties through interpolation from the fluid particles of neighboring subdomains. Specifically, at the interface between $Γ_{i}$ and $Γ_{i + 1}$ , buffer particles belonging to $Γ_{i}$ interpolate their physical quantities from the fluid particles in $Γ_{i + 1}$ , while buffer particles in $Γ_{i + 1}$ obtain their properties from fluid particles in $Γ_{i}$ . To minimize kernel truncation errors and ensure smooth data transfer across subdomain interfaces, the second-order interpolation technique of [75] is adopted. The interpolated physical quantities are subsequently used to update particle positions using a first-order Euler time-integration scheme. Particles are allowed to dynamically switch roles depending on their spatial location: (a) buffer particles become fluid particles when they move into a fluid subdomain, (b) fluid particles convert into buffer particles upon entering a buffer zone, and (c) particles leaving the buffer region are removed. To maintain a uniform particle distribution within buffer zones, particle shifting is applied following the procedure of [60], where concentration gradients are evaluated using nearby fluid particles together with auxiliary particles positioned on a uniform stencil beyond the computational domain.

To permit mass exchange between neighboring subdomains, the outer boundary of each buffer region is discretized using a Cartesian grid with a spacing of $Δ s_{i} = d p_{i}$ , defining a set of mass accumulation nodes A. The accumulated mass $m_{a}$ at each node increases incrementally through ${\dot{m}}_{a}$ until it exceeds the critical threshold ${(Δ s_{i})}^{3} / ρ_{0}$ . Once this condition is met, a new particle is inserted into the buffer region at a distance of $Δ s_{i}$ from the outer boundary. The instantaneous mass flux through each node is calculated as follows:

(11) ${\dot{m}}_{A} = \{\begin{matrix} max (0, - ρ_{m_{A}} v_{m_{A}} \cdot \bar{n} Δ s_{i}^{2} Δ t_{S P H}) & if {\dot{m}}_{A} < 0, \\ - ρ_{m_{A}} v_{m_{A}} \cdot \bar{n} Δ s_{i}^{2} Δ t_{S P H} & if {\dot{m}}_{A} \geq 0 . \end{matrix}$

where

ρ_{m_{A}}

and

v_{m_{A}}

are the interpolated density and velocity at node A,

\bar{n}

is the unit normal vector to the interface between adjacent subdomains, and

Δ t_{S P H}

is the SPH numerical timestep.

2.2. MBD Structural Solver

2.2.1. Governing Equations

Rigid body motions of fluid-driven floating structures, subjected to both mechanical constraints and mooring line forces, are computed using the multibody dynamics (MBD) module of Project Chrono [54]. Mechanical constraints used in AEM are resolved internally within Project Chrono, whereas moorings and flexible connectors are managed by MoorDyn+, as detailed in Section 2.3. The generalized positions $q$ of the floating bodies, driven by their velocities $v$ and the total applied forces $f_{t}$ , dynamically evolve in time according to:

(12) $\frac{d q}{d t} = L (q) v$

(13) $M \frac{d v}{d t} = f_{t} (t, q, v)$

(14) $f_{t} (t, q, v) = f_{e} - f_{c}$

where

L (q)

denotes the transformation matrix,

M

the mass matrix, and

f_{e}

and

f_{c}

the external and constraint force vectors, respectively.

The magnitude of the constraint force $f_{c}$ acting on each linear spring at a given timestep is expressed as:

(15) $∥ f_{c} ∥ = k (l - l_{0})$

where k is the spring stiffness, l is the instantaneous spring length, and

l_{0}

is its rest length. The components of

f_{c}

are evaluated relative to the attachment points of each spring on the floating bodies.

2.2.2. SPH-MBD Coupling

The coupling between DualSPHysics and Project Chrono follows the framework introduced by [55]. Within this coupling scheme, fluid-driven floating bodies are represented by a discrete set of $N_{k}$ boundary particles k, which are uniformly distributed with a spacing equal to $d p$ . These bodies are treated as rigid, preserving the initial distribution of their boundary particles throughout the simulation. Each solid particle interacts with neighboring fluid particles through the same kernel function used for fluid–fluid particle interactions and experiences the accumulated influence of fluid pressures, fluid viscous forces, and gravitational acceleration. The force per unit mass $f_{k}$ acting on an individual solid particle is thus evaluated as follows:

(16) $f_{k} = \sum_{b = 1}^{N_{b}} [m_{b} (\frac{P_{k} + P_{b}}{ρ_{k} ρ_{b}}) + Π_{k}] \nabla_{k} W_{k b} + g$

The summation of all forces per unit mass experienced by the solid particles are translated, within the communication interface of the two solvers, to total hydrodynamic linear force $F$ and angular force $T$ vectors acting on each rigid floating body. These two quantities, defined by Equations (17) and (18), represent the external forcing vector $f_{e}$ in Equation (14), which is applied at the centre of mass $R_{0}$ of each body:

(17) $F = M \frac{d V}{d t} = \sum_{k = 1}^{N_{k}} m_{k} f_{k}$

(18) $T = I \frac{d Ω}{d t} = \sum_{k = 1}^{N_{k}} m_{k} f_{k} \times (r_{k} - R_{0})$

where M is the total body mass derived from the material density

ρ_{s}

V

and

Ω

are the linear and angular velocity vectors, respectively,

I

is the mass moment of the inertia matrix, and

r_{k}

is the position vector of the

k^{th}

solid particle.

The equations of motion for rigid bodies (Equations (12)–(14)) are then integrated within Project Chrono using the hydrodynamic loads calculated using Equations (17) and (18). The MBD solver advances with its own internal timestep $t_{CH}$ , performing multiple iterations until the time alignment condition $t_{CH} \geq Δ t_{SPH}$ with the SPH solver is achieved. Upon reaching synchronization, the updated kinematic state of each body, comprised of linear velocity $V$ , angular velocity $Ω$ , and centre of mass position $R_{0}$ vectors, is transmitted back to the SPH solver. The corresponding solid particle properties are then updated within the fluid domain. This bidirectional coupling process is repeated throughout the simulation until the prescribed simulation time is reached.

2.3. Lumped-Mass Mooring Solver

2.3.1. Governing Equations

MoorDyn+ [56] is an enhanced version of the original MoorDyn model [57], developed to simulate mooring line dynamics using a lumped-mass formulation. In this approach, both mooring lines and flexible connectors are modelled within Moordyn+, where each element is discretized into $N_{m}$ nodes m, modeled as point masses, and consequently divided into $N_{m} - 1$ segments that connect adjacent nodes. Each segment is represented by spring–damper elements that reproduce the line’s axial viscoelastic behaviour, defined by constant physical and material properties such as its length $l_{M}$ , diameter $d_{M}$ , density $ρ_{M}$ , Young’s modulus $E_{M}$ , and internal damping coefficient $C_{M}$ . Bending and torsional effects are neglected in this formulation.

Therefore, the internal forces in each segment result from the combination of axial elasticity and damping, along with the effects of weight and buoyancy, where the submerged weight and buoyant forces are treated jointly. A segment located between nodes m and $m + 1$ is denoted by the subscript $m + \frac{1}{2}$ , while the segment located between nodes $m - 1$ and m is denoted by the subscript $m - \frac{1}{2}$ . The submerged weight at node m (Equation (20)) is expressed as the product of the submerged weights of the adjacent segments (Equation (19)):

(19) $W_{m + \frac{1}{2}} = A_{M} l_{M} (ρ_{0} - ρ_{M}) g$

(20) $W_{m} = \frac{1}{2} (W_{m + \frac{1}{2}} + W_{m - \frac{1}{2}}) {\hat{e}}_{z}$

where cross-sectional area

A_{M} = π {d_{M}}^{2} / 4

and

{\hat{e}}_{z}

is an upward unit vector. The axial tension

N

and damping force

C

vectors in each segment depend on the instantaneous strain

ε

between the nodes of each segment and are expressed as:

(21) $N_{m + \frac{1}{2}} = E_{M} A_{M} ε_{m + \frac{1}{2}}$

(22) $C_{m + \frac{1}{2}} = C_{M} A_{M} {\dot{ε}}_{m + \frac{1}{2}}$

(23) $ε_{m + \frac{1}{2}} = (\frac{∥ r_{m + 1} - r_{m} ∥}{l_{M}} - 1)$

where

r

is the position vector of nodes, and

\dot{(\cdot)} = \partial (\cdot) / \partial t

External forces acting on the mooring lines are defined by hydrodynamic drag and added mass. Seabed reaction forces are neglected in this study, although they can be modelled within the MoorDyn+ model, as the investigated validation cases do not consider such phenomena (Section 3). The drag components are evaluated following the Morison equation, split into normal $D_{n m}$ and tangential $D_{t m}$ contributions as follows:

(24) $D_{n m} = \frac{1}{2} ρ_{0} C_{d n} l_{M} d_{M} ∥ ({\dot{r}}_{m} \cdot {\hat{q}}_{m}) {\hat{q}}_{m} - {\dot{r}}_{m} ∥ [({\dot{r}}_{m} \cdot {\hat{q}}_{m}) {\hat{q}}_{m} - {\dot{r}}_{m}]$

(25) $D_{t m} = \frac{1}{2} ρ_{0} C_{d t} l_{M} d_{M} ∥ - ({\dot{r}}_{m} \cdot {\hat{q}}_{m}) {\hat{q}}_{m} ∥ [- ({\dot{r}}_{m} \cdot {\hat{q}}_{m}) {\hat{q}}_{m}]$

where

C_{d n}

and

C_{d t}

are the drag coefficients in the normal and tangential directions, respectively, and

{\hat{q}}_{m}

is the local tangent vector at node m. The added mass contribution is described through the matrix

a_{m}

(26) $a_{m} = ρ_{0} A_{M} l_{M} [C_{a n} (I_{m} - {\hat{q}}_{m} {\hat{q}}_{m}^{T}) + C_{a t} ({\hat{q}}_{m} {\hat{q}}_{m}^{T})]$

where

C a n

and

C a t

are the added-mass coefficients in the transverse and tangential directions, respectively, and

I_{m}

is the identity matrix.

Combining all contributions, the equation of motion governing each node becomes:

(27) $\underset{Mass and Added Mass}{\underset{⏟}{(m_{m} + a_{m})}} {\ddot{r}}_{m} = \underset{Axial Tension}{\underset{⏟}{(N_{m + \frac{1}{2}} - N_{m - \frac{1}{2}})}} + \underset{Axial Damping}{\underset{⏟}{(C_{m + \frac{1}{2}} - C_{m - \frac{1}{2}})}} + \underset{Weight}{\underset{⏟}{W_{m}}} + \underset{Drag}{\underset{⏟}{D_{n m} + D_{t m}}}$

where

m_{m}

is the node mass matrix and

\ddot{(\cdot)} = \partial^{2} (\cdot) / \partial t^{2}

2.3.2. SPH-MoorDyn+ Coupling

The coupling between DualSPHysics and MoorDyn+ follows the framework introduced by [52]. Within this scheme, SPH calculates the kinematic state of the floating body, namely its $V$ , $Ω$ and $R_{0}$ vectors, and transmits it to MoorDyn+ together with $Δ t_{S P H}$ . During each SPH timestep, MoorDyn+ employs a fixed internal timestep of $Δ t_{M} = 10^{- 5} s$ , which increments until $t_{M} \geq Δ t_{SPH}$ , thereby ensuring time alignment between the two solvers. Upon reaching synchronization, the updated kinematic state of each body, comprised of linear $d V / d t$ and angular $d Ω / d t$ acceleration vectors, is transmitted back to the SPH solver. The corresponding solid particle properties are then updated within the fluid domain. This bidirectional coupling process is repeated throughout the simulation until the prescribed simulation time is reached.

2.4. Applied Element Method

The AEM framework, initially developed for beam structures within the SPH–MBD coupling approach [40] and later extended to thin plates [41], is herein employed for the analysis of VFFS. In this formulation, a plate with dimensions $L_{x}$ , $L_{y}$ and H along the X-, Y-, and Z-axes, respectively, is discretized into a uniform grid of box-shaped rigid bodies (Figure 1). The origin of the coordinate system, for the cases defined by the AEM, is located at the midpoint of the upstream face of the floating structure relative to the wave direction. The discretization results in $N_{x} \times N_{y}$ rigid bodies, where $N_{x}$ and $N_{y}$ represent the number of divisions along the X- and Y-directions. Each body has dimensions $Δ_{x}$ , $Δ_{y}$ and H, given as follows:

(28) $Δ_{x} = \frac{L_{x}}{N_{x}}, Δ_{y} = \frac{L_{y}}{N_{y}} .$

The elastic interaction between adjacent rigid elements is modeled by sets of springs placed along each interface (Figure 1). Each interface contains sets of springs composed of one normal spring, which transfers axial forces, and two shear springs, which transfer tangential forces. The combined effects of these springs reproduce the stress field within a representative volume of the plate, highlighted with the orange colour in Figure 2 for a set of springs. Only in-plane normal stresses are considered, while through-thickness normal stresses are neglected. Discrete control points define the attachment positions of the spring sets along each interface. In the present formulation, the control points are located along the plate thickness and centered longitudinally on the interface (Figure 2), enabling a more accurate representation of through-thickness stress variation. The density of control points depends on the mesh resolution and on the dominant deformation mode of the elements, where shear-dominated behaviour or coarse meshes typically require a higher control-point density [76].

The stiffness of each spring, used in Equation (15), is determined from the material properties and the volume represented by each spring set:

(29) $k = \{\begin{matrix} \frac{E w d}{l}, & for normal springs \\ \frac{G w d}{l}, & for shear springs \end{matrix}$

where E is the Young’s modulus,

G = E / [2 (1 + ν)]

is the shear modulus, and

ν

is the Poisson’s ratio of the material. The geometric parameters l, w, and d denote the spring’s effective length, width, and representative thickness (Figure 2). For springs connecting elements along the X-axis,

l = Δ_{x}

and

w = Δ_{y}

, whereas for those along the Y-axis,

l = Δ_{y}

and

w = Δ_{x}

. The thickness d is user-defined and set to

d p

by default. For springs located along fixed boundaries, the corresponding l and w values are halved, while at free boundaries, they are increased by one half.

In the classical AEM formulation [77], the global stiffness matrix is assembled by summing the stiffness contributions of all springs, assuming fixed spring-force directions. This results in a constant stiffness matrix throughout the simulation, suitable for small-deformation analyses. To account for large non-linear displacements, [78] introduced a dynamic stiffness matrix formulation by geometrically updating the spring orientations at each time step. Consequently, the global stiffness matrix is continuously reformulated according to the evolving spring configurations. The present AEM implementation adopts this dynamic treatment of spring directionality, enabling accurate modelling of large-displacement and geometrically non-linear structural responses.

3. Validation

For such multidisciplinary applications, analytical solutions are not available, and the model’s validity is assessed using physical experiments. The coupled numerical model is therefore validated both quantitatively and qualitatively using two distinct experimental configurations: (a) an array of interconnected floating rigid pontoons (Section 3.1) and (b) a uniform flexible floating structure (Section 3.2). For each configuration, a detailed description of the experimental setup is presented (Section 3.1.1 and Section 3.2.1), followed by a description of its numerical formulation in the utilized numerical model (Section 3.1.2 and Section 3.2.2). The hydrodynamic response of the floating bodies obtained from the numerical simulations are then presented and validated against the experimental datasets in Section 3.1.3 and Section 3.2.3.

3.1. Moored Array of Interconnected Floating Rigid Pontoons

3.1.1. Experimental Setup

An array of six floating rigid pontoons, interconnected by linear soft connectors, was investigated under wave attack in [7]. The array was positioned at the centre of a wave basin with a length of $100 m$ , width of $3.8 m$ , and depth of $2.245 m$ , both longitudinally and transversally, to minimize wave reflections from the basin boundaries. A series of wave conditions were generated using a flap-type wavemaker, with the present study focusing on a regular wave case of amplitude $A_{w} = 0.0283 m$ and period $T_{w} = 1.37 s$ . Wave absorption was achieved using a dissipative beach profile with a length of $8.5 m$ , located at the opposite side of the wave basin relative to the wavemaker.

Each individual floating body had a length of $78.3 \times 10^{- 3} m$ , a width of $48.3 \times 10^{- 3} m$ , a thickness of $10 \times 10^{- 3} m$ , a dry material density of $313 {kg / m}^{3}$ , a wet material density of $352 {kg / m}^{3}$ , and a free-floating draft of $3.57 \times 10^{- 3} m$ . The bodies were arranged in a $3 \times 2$ configuration, with three of the pontoons positioned at an equal distance from the wavemaker along the longitudinal axis of the basin, while the remaining three were positioned with the same equal distance from the wavemaker, but further downstream. All pontoons maintained an equal spacing of $16.7 \times 10^{- 3} m$ in both the longitudinal and transverse directions, such that the geometrical centre of the array coincided with the geometrical centre of the basin. The outer corner bodies of the array were moored by four horizontal mooring lines located above the Still Water Level (SWL) to restrain the drift of the array. In the horizontal plane, each mooring line formed an angle of $135^{\circ}$ with the adjacent sides of the connected floating body, while the lines on the opposite side were attached to the lateral boundaries of the physical basin. In addition, all bodies were interconnected using soft connectors that shared the same geometric and material characteristics as the mooring lines, arranged in a mixed configuration (the reader is referred to [7] for a detailed illustration of the mooring and connector configurations). Both the mooring lines and the connectors consisted of lightweight wires with a material density of $5 \times 10^{- 4} kg / m$ and a diameter of $9 \times 10^{- 4} m$ . Furthermore, a pretension of $0.245 N$ was applied solely to the mooring lines.

Tensile forces in the mooring lines were measured using force sensors with a rated capacity of $200 N$ and a ratio error of $0.2 %$ , corresponding to an acceptable error of $0.4 N$ . Each sensor was connected to the mooring lines through a linear spring, which in order to represent the same restoring behaviour with the mooring system, had a stiffness of $k = 0.95 N / m$ . Model motions were captured by an optical tracking system composed of four cameras, tracking six retro-reflective markers mounted on lightweight pillars positioned at the centre of each float to minimize hydrodynamic interference. Calibration of the motion system ensured a measurement accuracy of $0.5 mm$ .

3.1.2. Numerical Setup

A digital twin of the experimental wave basin was developed to replicate the physical setup within a reduced and computationally optimized numerical domain. The basin dimensions, namely length, width, and depth, were proportionally scaled down to improve computational efficiency without compromising the repetitiveness of the physical experiments. Along the X-axis, the upstream boundary of the computational domain consists of a piston-type wavemaker positioned half a wavelength upstream of the moored array. Downstream of the moored array, at a distance of half a wavelength, a numerical damping zone was implemented to absorb incident waves. This damping region spans half a wavelength and is followed by a vertical wall that defines the downstream boundary of the computational domain. Along the Y-axis, lateral boundaries are also confined by vertical walls, defining the basin width equal to the array width extended by half a wavelength. The water depth was reduced to $1.67 m$ to further limit computational costs while maintaining deep-water conditions, according to the Le Méhauté graph [79]. According to this graph, deep-water conditions can be assumed when the ratio of water depth to wavelength exceeds the value of 0.5, indicating that the orbital motion of fluid particles in not distorted due to interactions with the basin’s bottom. Consequently, the modified water depth was selected to satisfy this criterion, ensuring that wave propagation remains unaffected, while also decreasing computational cost.

To further improve computational efficiency, the numerical domain was divided into four subdomains (Section 2.1.2). The first three, denoted as $Γ_{1}$ , $Γ_{2}$ , and $Γ_{3}$ , span the entire computational domain along the X (parallel to wave direction) and Y (in-plane perpendicular to wave direction) directions, but each represents a distinct portion of the water column along the Z-axis (out of plane perpendicular to wave direction). Herein, computational domain specifically refers to the region occupied by fluid or boundary particles, distinguishing it from the overall simulation domain. Subdomain $Γ_{1}$ corresponds to the deepest zone and employs the coarsest inter-particle resolution, whereas $Γ_{3}$ represents the free-surface zone. The resolution refinement follows the relation $d p_{3} = d p_{2} / 2 = d p_{1} / 4 \approx A_{w} / 3$ . The fourth subdomain $Γ_{4}$ encompasses the floating array region, maintaining a clearance around the structure in all directions, equal to the array’s length along the X-axis on each side, half the array’s width along the Y-axis on each side and half the height of subdomain $Γ_{3}$ along the Z-axis. This subdomain employs an inter-particle resolution of $d p_{4} = d p_{3} / 2 = H / 2$ . A schematic representation of the domain decomposition is provided in Figure 3, where subdomains $Γ_{1}$ – $Γ_{3}$ are shown with half their actual width and truncated in both the X and Z directions for illustration purposes. Only their bottom surfaces are depicted as planar, while subdomain $Γ_{4}$ is shown at its actual dimensions within the rectangular green box. The total simulation time was also limited to $t_{max} = 12 T_{w}$ .

Since the mooring lines are positioned above the SWL, the moorings are modeled without considering external hydrodynamic forces (Section 2.3.1). The geometric and material properties are consistent with the experimental setup. Both mooring lines and connectors are modelled using the mooring line equations (Section 2.3.1), where $N_{m} = 1$ and axial stiffness $E_{M} A_{M}$ of each segment is defined as $0.95 l_{M}$ , aligning with the spring stiffness used to connect the mooring line to the force sensors. As internal damping parameters are not provided from the experimental dataset, $C_{M} A_{M}$ is herein defined using a damping ratio of $1 %$ , solely to damp low-frequency vibrations of mooring segments and ensure numerical stability. The pretension of the mooring lines was also tuned to match the experimental model results, based on the mean values of the measured mooring tensions, resulting in a pretension of $0.265 N$ for the upstream mooring lines and $0.259 N$ for the downstream mooring lines. Pretension was applied by modifying the unstretched length of the mooring line $l_{M 0}$ , adjusted as a function of current line length $l_{M 1}$ , target pretension $N_{T}$ , and mooring line axial stiffness $E_{M} A_{M}$ , according to $l_{M 0} = l_{M 1} / (1 + (N_{T} / E_{M} A_{M}))$ .

3.1.3. Results

Prior to subjecting the array to regular wave excitation, an experimental heave decay test for a single pontoon was numerically reproduced without any moorings or connectors, with the objective of verifying the damping behaviour associated with fluid viscous forces (Equation (8)). In the experimental campaign, the pontoon was initially displaced by $- 1.3 \times 10^{- 3} m$ along the Z-axis, which led to a first heave peak of $3.1 \times 10^{- 3} m$ . This response is considered nonphysical, as the oscillation amplitude increases rather than decays. Accordingly, comparison with the experimental dataset is performed only after $t_{E X P} = 0.08 s$ , at which point the heave displacement is $3.1 \times 10^{- 3} m$ and is given as the initial displacement for the numerical heave decay test. As a result, Figure 4a depicts the heave response of a single pontoon as measured experimentally from $t_{E X P} = 0.08 s$ to $t_{E X P} = 0.50 s$ , and predicted numerically from $t_{S P H} = 0.00 s$ to $t_{S P H} = 0.42 s$ , shown with time-tick intervals of $0.10 s$ . An error of $9.23 %$ in the natural oscillation period was obtained, which was acceptable for investigating the heave response of the fully moored array, as repetitions in the experimental heave decay tests exhibit even greater variance.

Subsequently, an experimental surge decay test for the full moored array with all pontoons, moorings, and connectors was numerically reproduced to verify the restoring behaviour of the moored array. In this simulation, all pontoons were displaced by $3.4 \times 10^{- 1} m$ along the X-axis. As a result, Figure 4b depicts the surge response of a single pontoon within the moored array, as measured experimentally from $t_{E X P} = 0.41 s$ (corresponding to the starting time frame of the surge decay test) to $t_{E X P} = 4.41 s$ and predicted numerically from $t_{S P H} = 0.00 s$ to $t_{S P H} = 4.00 s$ , shown with time-tick intervals of $1.00 s$ . The moored array showcased an underdamped numerical response relative to the experimental response. Other numerical investigations [16,18] replicating the same experimental configuration consistently report similar underdamped numerical responses. This discrepancy is primarily attributed to an inadequate representation of the mooring lines’ damping characteristics, which are herein almost neglected due to lack of experimental documentation on their material properties. Additional deviations may also be attributed to errors inherent in experimental measurements. In previous numerical studies, tuning of viscous fluid force coefficients has been employed to enable reliable assessment of surge motions under wave excitation. In the present work, such parameter tuning is intentionally omitted; therefore, a comparison of surge motions under wave excitation is not presented.

Prior to evaluating the hydrodynamic response of the moored array, a sensitivity analysis of particle refinement ratios between adjacent subdomains was conducted in an empty basin, for $d p_{i} = 1.5 d p_{i + 1}$ , $d p_{i} = 2.0 d p_{i + 1}$ and $d p_{i} = 2.5 d p_{i + 1}$ , where $i = 1, . . ., 4$ and $d p_{4} = H / 2$ . Convergence was assessed in Figure 5 by comparing fluid surface elevations calculated analytically from 2nd order Stokes wave theory and predicted numerically from $t_{S P H} = 3.39 s$ to $t_{S P H} = 10.34 s$ . Values are shown with time-tick intervals of $T_{w} = 1.37 s$ . Wave amplitude error values for $d p_{i} = 1.5 d p_{i + 1}$ , $d p_{i} = 2.0 d p_{i + 1}$ , and $d p_{i} = 2.5 d p_{i + 1}$ were $2.06 %$ , $2.58 %$ , and $4.66 %$ , respectively, while wave period error values were consistently below $1.00 %$ for all cases. Thus, the selected refinement of $d p_{i} = 2.0 d p_{i + 1}$ was confirmed to showcase great prediction accuracy within a feasible computational cost. This analysis simultaneously verified the accuracy of wave generation and absorption in the SPH model.

The heaving motion of the entire moored array is given for a representative pontoon within the array, as only minor variations were observed among individual units. As a result, Figure 6a depicts the heave response of the representative pontoon, as measured experimentally from $t_{E X P} = 71.00 s$ to $t_{E X P} = 77.46 s$ and predicted numerically from $t_{S P H} = 3.39 s$ to $t_{S P H} = 9.85 s$ , shown at time-tick intervals of $T_{w} = 1.37 s$ . The comparison indicates a mean oscillation amplitude error of $2.63 %$ and an oscillation period error below $1 %$ . Furthermore, the mooring line tensions exhibited nearly identical maximum and minimum values (Figure 6b) compared with the experimental results for both the upstream and downstream mooring lines.

The variable resolution scheme utilized in this validation case required 4,519,660 fluid particles and 27.67 GPU-hours on an NVIDIA GeForce RTX 3090 GPU for a simulation duration of $12 T_{w}$ . In contrast, a uniform resolution scheme would require 57,328,128 fluid particles, with the SPH model predicting an estimated 255.80 GPU-hours for the same simulation duration. This corresponds to a 92.1% reduction in fluid particles and an estimated 89.2% reduction in computational time achieved by the variable resolution scheme.

3.2. Moored Very Flexible Floating Structure

3.2.1. Experimental Setup

A continuous moored very flexible floating structure was investigated under regular wave attack in [80], at a wave basin with a length of $80 m$ , width of $2.75 m$ , and water depth of $1.00 m$ . Two regular wave conditions were generated using a flap-type wavemaker, with the present study focusing on the regular wave case of amplitude $A_{w} = 0.02 m$ and period $T_{w} = 0.796 s$ . Wave absorption was achieved using a dissipative beach profile, located at the opposite side of the wave basin relative to the wavemaker.

The floating structure consisted of a closed-pore neoprene foam rubber sheet with a length of $L_{x} = 4.95 m$ , width of $L_{y} = 1.02 m$ , and thickness of $H = 5 \times 10^{- 3} m$ . The material properties were characterized by a density of $ρ_{s} = 116 {kg / m}^{3}$ , a Young’s modulus of $E = 5.60 \times 10^{5} Pa$ , and a Poisson’s ratio of $ν = 0.4$ . All corners of the floating structure were moored by four horizontal mooring lines located at the SWL to restrain its drift. In the horizontal plane, each mooring line formed an angle of $100 . 08^{\circ}$ with the short adjacent side of the floating structure and an angle of $169 . 92^{\circ}$ with the long adjacent side, while the lines on the opposite side were attached to the lateral boundaries of the physical basin. The upstream mooring lines, towards the wavemaker, were stiff DYNEEMA mooring lines with a diameter of $0.2 \times 10^{- 3} m$ and installed without pretension. The downstream mooring lines, towards the dissipative beach, were soft lines with a diameter of less than $1 \times 10^{- 3} m$ . Furthermore, a pretension of $0.1 N$ was applied solely to the downstream mooring lines, with an axial stiffness $E_{M} A_{M} = 0.475 N$ at the pretensioned elongation.

Displacements of the floating structure under wave attack were determined using a digital image correlation system, where images of the undeformed and deformed configurations were compared to derive surface displacements. A two-camera setup with a resolution of 4 megapixels was used, in which one pixel corresponded to $7 \times 10^{- 4} m$ . Correlation between the deformed and undeformed images was achieved through a random speckle pattern applied to the structure surface, with an average speckle size of $9 \times 10^{- 3} m$ . Calibration was carried out with a $1 m \times 1 m$ target with $7.5 \times 10^{- 3} m$ grid spacing, ensuring accurate transformation of the images into the real-world coordinate system for the reconstruction of surface displacements.

3.2.2. Numerical Setup

A digital twin of the experimental wave basin was developed following a formulation similar to that described in Section 3.1.2. Along the X-axis, the upstream boundary of the computational domain features a piston-type wavemaker located one wavelength upstream of the flexible floating structure. Downstream of the structure, at a distance of one wavelength, a numerical damping zone of one wavelength in length was implemented, followed by a vertical wall defining the downstream boundary. Along the Y-axis, the lateral boundaries are confined by vertical walls, giving the basin a width equal to that of the floating structure extended by approximately half a wavelength. To reduce computational cost while maintaining deep-water conditions [79], the water depth was set to $0.55 m$ (as described in Section 3.1.2).

The four numerical subdomains were constructed using the same formulation strategy described in Section 3.1.2. The mesh refinement followed the relations $d p_{3} = d p_{2} / 2 = d p_{1} / 4 = A_{w} / 2$ and $d p_{4} = d p_{3} / 2 = H / 2$ . For Subdomain $Γ_{4}$ , a clearance around the structure in all directions was maintained, equal to $0.20 m$ along the X-axis on each side, $0.10 m$ along the Y-axis on each side, and $0.15 m$ along the Z-axis. A schematic representation of the domain decomposition is provided in Figure 7, where Subdomains $Γ_{1}$ – $Γ_{3}$ are shown with half their actual width and truncated in both the X and Z directions for illustration purposes. Only their bottom surfaces are depicted as planar, while Subdomain $Γ_{4}$ is shown at its actual dimensions within the rectangular green box. The total simulation time was also limited to $t_{max} = 20 T_{w}$ .

In this case, since the mooring lines are positioned at the SWL, the moorings are modeled considering external hydrodynamic forces (Section 2.3.1) using $C_{d n} = 1.20$ , $C_{d t} = 0.20$ , $C_{a n} = 1.00$ and $C_{a t} = 0$ , with internal damping corresponding to $1 %$ of critical damping, used solely to damp low-frequency vibrations of mooring segments and ensure numerical stability, and $N_{m} = 1$ . The geometric and material properties match those of the experimental setup (with pretension applied as described in Section 3.1.2), while the axial stiffness of the upstream mooring lines was set to $E_{M} A_{M} = 66.71 N$ , consistent with DYNEEMA products of similar geometric characteristics [81].

3.2.3. Results

Numerical vertical displacements $d_{Z}$ , along the Z-axis, normalized by the incident wave amplitude $A_{w}$ , were compared qualitatively and quantitatively with the experimental measurements. The full top-surface elevation fields of the floating structure along the X- and Y-axes are presented in Figure 8. To mitigate edge-related experimental measurement fluctuations, the data were cropped by $0.5 m$ on each side along the Y-axis and $0.5 m$ from the downstream edge along the X-axis. The comparison shows that the numerical model accurately reproduces the spatial distribution of crests and troughs in the surface displacements and captures the relative amplitude variations across the structure. To further validate the model, three cross-sections along the Y-axis, namely at $Y = - 0.4 m$ , $Y = 0 m$ , and $Y = 0.4 m$ (indicated by red dashed lines in Figure 8), are plotted in Figure 8. The results show excellent agreement at the initial peaks, with the numerical solution exhibiting a mild decay near the downstream boundary.

Quantitatively, vertical displacements $d_{Z}$ normalized by the incident wave amplitude $A_{w}$ confirm a strong agreement between the numerical and experimental results (Figure 9a), with RMSE values of $8.22 %$ , $8.96 %$ , and $8.88 %$ for $Y = - 0.4 m$ , $Y = 0.0 m$ , and $Y = 0.4 m$ , respectively. The oscillation amplitude errors per wave period further indicate consistent model performance, remaining within $\pm 10 %$ across all sections. The slightly negative amplitude errors observed at $Y = - 0.4 m$ and $Y = 0.4 m$ suggest a minor underestimation of oscillation intensity near the structure edges, while the near-zero deviations at $Y = 0.0 m$ confirm accurate prediction in the central region. For the downstream oscillation amplitudes, a greater underestimation is observed, attributed to inherent energy dissipation in SPH during wave propagation.

To quantify this dissipative behaviour of propagating waves, the simulation was repeated in an empty basin, where the flexible floating structure was removed from the simulation domain. The normalized fluid surface elevations $η$ , over the target wave amplitude $A_{w}$ , are shown in Figure 9b, recorded along the X-axis at the location where the structure was previously situated and along the Y-axis at $Y = 0.0 m$ . An average peak reduction value of $4.45 %$ was observed over four wavelengths, with the maximum reduction occurring at the last peak with a value of $8.07 %$ . Overall, the results validate the model’s ability to capture both the spatial and temporal characteristics of the wave-induced displacements with high fidelity.

The variable resolution scheme utilized in this validation case required 11,587,708 fluid particles and 38.61 GPU-hours on an NVIDIA GeForce RTX 3090 GPU, for a simulation duration of $20 T_{w}$ . In contrast, a uniform resolution scheme would require 52,968,960 fluid particles, with the SPH model predicting an estimated 139.29 GPU-hours for the same simulation duration. This corresponds to a 78.1% reduction in fluid particles and an estimated 72.3% reduction in computational time achieved by the variable resolution scheme.

4. Discussion

The proposed three-dimensional SPH framework was successfully validated for accurately predicting the hydrodynamic performance of moored modularized floating structures and VFFS. The initial validation case focused on modularized structures, assessing the combined use of the SPH–MBD and SPH-MoorDyn+ coupling frameworks for complex multibody systems, including floating pontoons, their coupled hydrodynamic behaviour through flexible connectors, and mooring lines. Although the proposed SPH model does not require parameter tuning, its performance is sensitive to the artificial viscosity parameter (Equation (8)), which is empirically set to a default value, showcasing an acceptable fluid–structure interaction behaviour for the current test case (Figure 4a). However, a more accurate physical dissipation formulation in the SPH equations could further enhance the accuracy of the results. For regular wave conditions, heave (Figure 6a) response was validated with high accuracy, indicating that the governing equations were correctly formulated and that no further wave-specific parameter tuning was necessary. Furthermore, mooring line tensions were accurately predicted (Figure 6b), demonstrating the robustness of the SPH–MoorDyn+ coupling, even when variable resolution schemes were employed. Yet, a major limitation of the proposed SPH approach was highlighted in the surge decay test (Figure 4b), where an accurate definition of damping characteristics for mooring lines is required.

A secondary validation case focused on moored VFFS structures, further confirming the accuracy of the AEM scheme within both the SPH-MBD and the SPH-MoorDyn+ coupling frameworks, as well as in their combined application under a variable resolution scheme. The numerical results demonstrated consistent agreement with experimental measurements in both the longitudinal and transverse directions of the structure (Figure 8), with significant discrepancies appearing only towards the downstream end of the structure (Figure 9a). These deviations are attributed to wave decay effects inherent to SPH-based wave propagation models, where non-physical energy dissipation leads to a gradual reduction in wave amplitude over several wavelengths (Figure 9b). This non-physical energy dissipation has a measurable impact both on long-wave propagation and the wave-induced hydroelastic response of VFFS, thereby limiting the applicability of the proposed SPH approach to simulation domains spanning only a few wavelengths.

For both validation cases, a variable resolution scheme was utilized, conducting simulations that would otherwise be computationally prohibitive, both in terms of runtime and data storage. By locally refining particle resolution only within targeted regions of the simulation domain, the SPH approach preserves high spatial fidelity only in the areas of interest for each specific case-study.

5. Conclusions

A versatile three-dimensional numerical approach has been developed and validated for the simulation of modularized and flexible floating structures subjected to wave loading. The proposed model integrates four key components, namely non-linear fluid dynamics, multibody dynamics, structural flexibility, and mooring line dynamics, within a unified computational environment. The combined utilization of SPH-MBD coupling, the Applied Element Method (AEM), and the SPH-MoorDyn+ mooring model enables accurate and efficient simulation of complex fluid–structure–mooring interactions.

The validation studies, involving both moored arrays of rigid interconnected pontoons and moored highly flexible floating structures, demonstrated that the framework can reliably predict hydrodynamic responses, structural deformations, and mooring line tensions under regular wave conditions. The coupling was shown to capture non-linear hydroelastic responses without requiring case-specific parameter tuning, relying solely on generally adapted default formulations for artificial viscosity and density diffusion, highlighting the model’s robustness and generality. However, the results also revealed two primary discrepancies in the proposed SPH method, its reliance on an empirically tuned artificial viscosity coefficient and an unphysical energy dissipation in propagating waves. Regarding the validation cases, the coupled floating-mooring system was not validated in any of the surge or sways degrees of freedom, while for fluid–structure interactions, only weakly non-linear phenomena were considered.

Future work should therefore focus on: (a) developing gradient kernel correction schemes that might limit non-physical energy dissipation in propagating waves [82]; (b) demonstrating the reliability of the proposed method in predicting horizontal mooring restoring forces; (c) investigating scenarios involving violent fluid flows, where DualSPHysics excels, interacting with geometrically non-linear floating structures; (d) implementing non-linear modelling of mooring dynamics; and (e) integration of hyperelastic material formulations within the AEM framework.

Overall, the developed SPH-based multiphysics framework provides a powerful numerical tool for the analysis and design of emerging floating modular and flexible systems in fully three-dimensional marine environments through the use of a variable resolution scheme. Its versatility allows extension to more complex scenarios without any parameter tuning, including non-linear irregular sea states, hinged connectors, and shared mooring systems, supported by relevant validated functionalities available in the DualSPHysics, Project Chrono, and MoorDyn+ coupled models.

Author Contributions

Conceptualization, R.I., V.S., E.L., and P.T.; methodology, R.I., V.S., and E.L.; software, R.I.; validation, R.I.; formal analysis, R.I.; investigation, R.I.; resources, R.I.; data curation, R.I.; writing—original draft preparation, R.I.; writing—review and editing, V.S., E.L., and P.T.; visualization, R.I.; supervision, P.T.; project administration, P.T.; funding acquisition, V.S., E.L., and P.T. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The original contributions presented in this study are included in the article. Further inquiries can be directed to the corresponding author.

Acknowledgments

The computational resources used in this work were provided by the VSC (Vlaams Supercomputer Centrum), funded by FWO and the Flemish Government. Furthermore, during the preparation of this manuscript, the authors used OpenAI (https://openai.com) to correct grammatical errors and enhance readability in the text. The authors have reviewed and edited the output and take full responsibility for the content of this publication.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures

Figure 1 Schematic illustration of a thin plate partitioned into rigid bodies for $N_{x} = N_{y} = 5$ , where each resulting rigid body has dimensions $Δ_{x}$ × $Δ_{y}$ × H. The highlighted region with red corresponds to the detail shown in Figure 2.

Figure 2 Detailed view of the interface between adjacent rigid bodies, illustrating four sets of springs distributed along the plate thickness. Each set of springs is composed of one normal spring (shown in red colour) and two shear springs (shown in blue and green colours). The region highlighted in orange indicates the representative volume associated with one set of springs.

Figure 3 Schematic representation of a moored array of six interconnected rigid pontoons within an SPH variable resolution scheme. The illustrated subdomains $Γ_{1}$ – $Γ_{3}$ are halved in width, truncated in the X- and Z-axes, and depicted by their bottom surfaces shown as planar. Subdomain $Γ_{4}$ is depicted at its full scale highlighted by the rectangular green box.

Figure 4 Comparison between numerical and experimental results for: (a) heave motions; and (b) surge motions of a free-floating individual pontoon. The experimental values are shown as black dashed lines, while the numerical values are shown as orange solid lines.

Figure 5 Comparison between numerical and analytical fluid surface elevations $η$ in an empty basin. Numerical results are shown for different particle refinement ratios between adjacent subdomains, namely $d p_{i} = 1.5 d p_{i + 1}$ , $d p_{i} = 2.0 d p_{i + 1}$ , and $d p_{i} = 2.5 d p_{i + 1}$ , where $i = 1, . . ., 4$ and $d p_{4} = H / 2$ . The analytical values are shown as black dashed lines, while the numerical values are shown as coloured solid lines.

Figure 6 Comparison between numerical and experimental results for: (a) Heave motions of an individual pontoon within a moored array of six interconnected rigid pontoons. The experimental values are shown as black dashed lines, while the numerical values are shown as orange solid lines. (b) Mooring line tensions for a moored array of six interconnected rigid pontoons for two representative mooring lines located upstream and downstream. The experimental responses are shown as black solid lines, while the numerical values are shown as orange solid lines.

Figure 7 Schematic representation of a moored very flexible floating structure within an SPH variable resolution scheme. The illustrated Subdomains $Γ_{1}$ – $Γ_{3}$ are halved in width, truncated in the X- and Z- axes, and depicted by their bottom surfaces shown as planar. Subdomain $Γ_{4}$ is depicted at its full scale highlighted by the rectangular green box.

Figure 8 Comparison between experimental (left) and numerical (right) normalized vertical displacements $d_{Z}$ , over the incident wave amplitude $A_{w}$ , on the top surface of the moored flexible floating structure. Structural coordinates are provided from 0.0 m to 4.95 m along the X-axis and from $- 0.45$ m to $0.45$ m along the Y-axis.

Figure 9 (a) Comparison between experimental and numerical normalized vertical displacements $d_{Z}$ , over the incident wave amplitude $A_{w}$ , on the top surface of the moored flexible floating structure. Structural coordinates are provided from 0.05 m to 4.8 m along the X-axis and at three transverse planes located at $Y = - 0.4$ m, $Y = 0.0$ m, and $Y = 0.4$ m. The experimental values are shown as black dashed lines, while the numerical values are shown as coloured dashed lines. (b) Numerical normalized fluid surface elevations $η$ , over the incident wave amplitude $A_{w}$ , in an empty wave basin without the moored flexible floating structure. Surface elevations are provided from 0.985 m to 5.985 m along the X-axis, and at the transverse plane located at $Y = 0.0$ m. The numerical values are shown as a black solid line, while the red dashed lines indicate the target wave envelope.

References

1. Claus, R.; López, M. Key issues in the design of floating photovoltaic structures for the marine environment. Renew. Sustain. Energy Rev.; 2022; 164, 112502. [DOI: https://dx.doi.org/10.1016/j.rser.2022.112502]

2. Wang, J.; Lund, P.D. Review of recent offshore photovoltaics development. Energies; 2022; 15, 7462. [DOI: https://dx.doi.org/10.3390/en15207462]

3. López, M.; Claus, R.; Soto, F.; Hernández-Garrastacho, Z.A.; Cebada-Relea, A.; Simancas, O. Advancing offshore solar energy generation: The HelioSea concept. Appl. Energy; 2024; 359, 122710. [DOI: https://dx.doi.org/10.1016/j.apenergy.2024.122710]

4. HelioRec. Wind and Wave-Resistant Floating Solar Technology. Available online: https://www.heliorec.com/tech (accessed on 1 October 2025).

5. Moss Maritime. Offshore Floating Solar Power. Available online: https://www.mossww.com/offshore-floating-solar-power/ (accessed on 1 October 2025).

6. Oceans of Energy. North Sea 3. Available online: https://oceansofenergy.blue/north-sea-3/ (accessed on 1 October 2025).

7. Jiang, Z.; Dai, J.; Saettone, S.; Tørå, G.; He, Z.; Bashir, M.; Souto-Iglesias, A. Design and model test of a soft-connected lattice-structured floating solar photovoltaic concept for harsh offshore conditions. Mar. Struct.; 2023; 90, 103426. [DOI: https://dx.doi.org/10.1016/j.marstruc.2023.103426]

8. Wang, B.; Lian, J.; Gao, X.; Liu, L.; Yao, Y.; Li, Y.; Meng, H. Coupled dynamic analysis of polyester flexible connector for a novel lattice-type floating photovoltaic array. Ocean Eng.; 2025; 320, 120354. [DOI: https://dx.doi.org/10.1016/j.oceaneng.2025.120354]

9. Lian, J.; Cao, Z.; Lu, W.; Li, P.; Liu, D.; Yao, Y.; Shao, N.; Li, H. Experimental investigation on the dynamic characteristics of connector and mooring force for a novel pontoon-type offshore floating photovoltaic array. Sol. Energy; 2025; 301, 113912. [DOI: https://dx.doi.org/10.1016/j.solener.2025.113912]

10. Li, Z.; Chen, D.; Feng, X. Hydroelastic and expansibility analysis of a modular floating photovoltaic system with multi-directional hinge connections. Ocean Eng.; 2023; 289, 116218. [DOI: https://dx.doi.org/10.1016/j.oceaneng.2023.116218]

11. Jiang, C.; Li, Q.; Hu, P.; Zhang, G.; Jiang, S.; el Moctar, O. Experimental analysis of wave–structure interaction for very large modular floating structures. Ocean Eng.; 2025; 315, 119763. [DOI: https://dx.doi.org/10.1016/j.oceaneng.2024.119763]

12. Mohapatra, S.C.; Amouzadrad, P.; Bispo, I.B.d.S.; Guedes Soares, C. Hydrodynamic response to current and wind on a large floating interconnected structure. J. Mar. Sci. Eng.; 2025; 13, 63. [DOI: https://dx.doi.org/10.3390/jmse13010063]

13. Ikhennicheu, M.; Blanc, A.; Danglade, B.; Gilloteaux, J.-C. OrcaFlex modelling of a multi-body floating solar island subjected to waves. Energies; 2022; 15, 9260. [DOI: https://dx.doi.org/10.3390/en15239260]

14. Song, J.; Imani, H.; Yue, J.; Yang, S. Hydrodynamic characteristics of floating photovoltaic systems under ocean loads. J. Mar. Sci. Eng.; 2023; 11, 1813. [DOI: https://dx.doi.org/10.3390/jmse11091813]

15. Zheng, Z.; Jin, P.; Huang, Q.; Zhou, B.; Xiang, R.; Zhou, Z.; Huang, L. Motion response and energy harvesting of multi-module floating photovoltaics in seas. Ocean Eng.; 2024; 310, 118760. [DOI: https://dx.doi.org/10.1016/j.oceaneng.2024.118760]

16. Baniya, J.; Dai, J.; Jiang, Z.; Zhang, C. Numerical study of soft-connected modular offshore floating photovoltaic array. Mar. Struct.; 2025; 104, 103893. [DOI: https://dx.doi.org/10.1016/j.marstruc.2025.103893]

17. He, W.; Liu, Z.; Yin, D.; Du, J.; Chang, A. Experimental and numerical research on the dynamic response of multi-module floating photovoltaics array. Ocean Eng.; 2025; 338, 121969. [DOI: https://dx.doi.org/10.1016/j.oceaneng.2025.121969]

18. Saitov, K.; Jiang, Z. Advancing offshore solar energy generation: Numerical modelling and calibration of a floating solar array against model tests. Sol. Energy; 2025; 288, 113225. [DOI: https://dx.doi.org/10.1016/j.solener.2024.113225]

19. Wang, B.; Li, Y.; Huang, L.; Yao, Y.; Qin, Y. Dynamic analysis of a novel star-type floating photovoltaic system with flexible connectors. Ocean Eng.; 2024; 304, 117854. [DOI: https://dx.doi.org/10.1016/j.oceaneng.2024.117854]

20. Lee, J.-H.; Paik, K.-J.; Lee, S.-H.; Hwangbo, J.; Ha, T.-H. Experimental and numerical study on the characteristics of motion and load for a floating solar power farm under regular waves. J. Mar. Sci. Eng.; 2022; 10, 565. [DOI: https://dx.doi.org/10.3390/jmse10050565]

21. Pan, Y.; Hong, J.; Xue, D.; Zhao, X.; Luo, M.; Liu, C.; Li, L. A SPH method to assess the dynamic behavior for multiconnected floating photovoltaic units consisting of non-uniform supporting frames in waves. Ocean Eng.; 2025; 323, 120622. [DOI: https://dx.doi.org/10.1016/j.oceaneng.2025.120622]

22. Li, X.; Xiao, Q.; Zhou, Y.; Ning, D.; Incecik, A.; Nicoll, R.; McDonald, A.; Campbell, D. Coupled CFD–MBD numerical modeling of a mechanically coupled WEC array. Ocean Eng.; 2022; 256, 111541. [DOI: https://dx.doi.org/10.1016/j.oceaneng.2022.111541]

23. Chen, H.; Almeida Medina, T.; Cercos-Pita, J.L. CFD simulation of multiple moored floating structures using OpenFOAM: An open-access mooring restraints library. Ocean Eng.; 2024; 303, 117697. [DOI: https://dx.doi.org/10.1016/j.oceaneng.2024.117697]

24. Tagliafierro, B.; Karimirad, M.; Altomare, C.; Göteman, M.; Martínez-Estévez, I.; Capasso, S.; Domínguez, J.M.; Viccione, G.; Gómez-Gesteira, M.; Crespo, A.J.C. Numerical validations and investigation of a semi-submersible floating offshore wind turbine platform interacting with ocean waves using an SPH framework. Appl. Ocean Res.; 2023; 141, 103757. [DOI: https://dx.doi.org/10.1016/j.apor.2023.103757]

25. Kong, M.; Zhang, X.; Ji, R.; Wu, H.; Yin, M.; Liu, H.; Sun, K.; Reabroy, R. Effects of wave–current interaction on hydrodynamic performance and motion response of a floating tidal stream turbine. J. Mar. Sci. Eng.; 2025; 13, 1520. [DOI: https://dx.doi.org/10.3390/jmse13081520]

26. Tang, Y.; Chen, X.; Huang, G.; Chen, M.; Zhang, K.; Jiang, Y.; Li, T.; Tao, T.; Yang, Y. Dynamic analysis of multi-module floating photovoltaic platforms with composite mooring system by considering tidal variation and platform configuration. Ocean Eng.; 2024; 312, 119243. [DOI: https://dx.doi.org/10.1016/j.oceaneng.2024.119243]

27. Zhang, G.; Bao, F.; Jiang, C.; el Moctar, O.; Xu, Z.; Zhou, B. Hydrodynamic performance of floating solar arrays with articulated modules for optimal wave adaptation. Appl. Ocean Res.; 2025; 164, 104762. [DOI: https://dx.doi.org/10.1016/j.apor.2025.104762]

28. Ricci, F.; Vacondio, R.; Tafuni, A. Multiscale smoothed particle hydrodynamics based on a domain-decomposition strategy. Comput. Methods Appl. Mech. Eng.; 2024; 418, 116500. [DOI: https://dx.doi.org/10.1016/j.cma.2023.116500]

29. Ricci, F.; Vacondio, R.; Domínguez, J.M.; Tafuni, A. Three-dimensional variable resolution for multi-scale modeling in smoothed particle hydrodynamics. Comput. Phys. Commun.; 2025; 313, 109609. [DOI: https://dx.doi.org/10.1016/j.cpc.2025.109609]

30. Ocean Sun. The Innovation. Available online: https://oceansun.no/the-innovation/ (accessed on 1 October 2025).

31. MARINBluewater Energy Services B.V.Endures B.V.Genap B.V.Jifmar Offshore ServicesTNO. Solar@Sea III Project Overview. Available online: https://www.marin.nl/en/jips/solarsea-iii (accessed on 1 October 2025).

32. Zhang, M.; Schreier, S. Review of wave interaction with continuous flexible floating structures. Ocean Eng.; 2022; 264, 112404. [DOI: https://dx.doi.org/10.1016/j.oceaneng.2022.112404]

33. Brown, S.A.; Xie, N.; Hann, M.R.; Greaves, D.M. Investigation of wave-driven hydroelastic interactions using numerical and physical modelling approaches. Appl. Ocean Res.; 2022; 129, 103363. [DOI: https://dx.doi.org/10.1016/j.apor.2022.103363]

34. Jiang, C.; el Moctar, O.; Zhang, G. Nonlinear modeling of wave–structure interaction for a flexible floating structure. Ocean Eng.; 2024; 300, 117489. [DOI: https://dx.doi.org/10.1016/j.oceaneng.2024.117489]

35. Jia, A.; Wang, X.; Li, A.; Liu, Y.; Jiang, S. Numerical investigation of the effects of edge constraints on the wave-induced hydroelastic responses of the single large floating module. Ocean Eng.; 2025; 326, 120885. [DOI: https://dx.doi.org/10.1016/j.oceaneng.2025.120885]

36. Jin, Z.; Fang, H.; Liu, Y. Numerical investigation on hydro-elastoplastic responses of floating laminated structures under wave action. Mar. Struct.; 2025; 104, 103884. [DOI: https://dx.doi.org/10.1016/j.marstruc.2025.103884]

37. Zhang, G.; Li, Q.; Jiang, C.; el Moctar, O.; Sun, Z. Constraint effects on the hydroelasticity of very large floating structures. Ocean Eng.; 2025; 331, 121304. [DOI: https://dx.doi.org/10.1016/j.oceaneng.2025.121304]

38. Ma, C.; Iijima, K.; Oka, M. Nonlinear waves in a floating thin elastic plate, predicted by a coupled SPH and FEM simulation and by an analytical solution. Ocean Eng.; 2020; 204, 107243. [DOI: https://dx.doi.org/10.1016/j.oceaneng.2020.107243]

39. Zhan, Y.; Luo, M.; Khayyer, A. An enhanced numerical wave tank for wave–structure interaction based on DualSPHysics+. Ocean Eng.; 2025; 340, 122413. [DOI: https://dx.doi.org/10.1016/j.oceaneng.2025.122413]

40. Ioannou, R.; Stratigaki, V.; Loukogeorgaki, E.; Troch, P. Establishing AEM structural framework within SPH-MBD coupling for hydro-viscoelastic response of very flexible floating structures. Appl. Ocean Res.; 2025; 164, 104760. [DOI: https://dx.doi.org/10.1016/j.apor.2025.104760]

41. Ioannou, R.; Stratigaki, V.; Loukogeorgaki, E.; Troch, P. Thin floating plate formulation in smoothed particle hydrodynamics: Extension to a three-dimensional applied element method framework. Coastal Eng.; 2025.submitted

42. Lyu, H.-G.; Sun, P.-N.; Huang, X.-T.; Zhong, S.-Y.; Peng, Y.-X.; Jiang, T.; Ji, C.-N. A review of SPH techniques for hydrodynamic simulations of ocean energy devices. Energies; 2022; 15, 502. [DOI: https://dx.doi.org/10.3390/en15020502]

43. O’Connor, J.; Rogers, B.D. A fluid–structure interaction model for free-surface flows and flexible structures using smoothed particle hydrodynamics on a GPU. J. Fluids Struct.; 2021; 104, 103312. [DOI: https://dx.doi.org/10.1016/j.jfluidstructs.2021.103312]

44. Shimizu, Y.; Khayyer, A.; Gotoh, H. An implicit SPH-based structure model for accurate fluid–structure interaction simulations with hourglass control scheme. Eur. J. Mech. B Fluids; 2022; 96, pp. 122-145. [DOI: https://dx.doi.org/10.1016/j.euromechflu.2022.07.007]

45. Khayyer, A.; Shimizu, Y.; Lee, C.H.; Gil, A.; Gotoh, H.; Bonet, J. An improved updated Lagrangian SPH method for structural modelling. Comput. Part. Mech.; 2024; 11, pp. 1055-1086. [DOI: https://dx.doi.org/10.1007/s40571-023-00673-z]

46. Bao, T.; Hu, J.; Wang, S.; Huang, C.; Yu, Y.; Shakibaeinia, A. An entirely SPH-based FSI solver and numerical investigations on hydrodynamic characteristics of the flexible structure with an ultra-thin characteristic. Comput. Methods Appl. Mech. Eng.; 2024; 431, 117255. [DOI: https://dx.doi.org/10.1016/j.cma.2024.117255]

47. Shi, Y.; Li, S.; Peng, S.; Chen, H.; He, M. Numerical modeling of flexible floating boom using a coupled SPH–FEM model. Coast. Eng. J.; 2018; 60, pp. 140-158. [DOI: https://dx.doi.org/10.1080/21664250.2018.1460060]

48. Hermange, C.; Oger, G.; Le Chenadec, Y.; Le Touzé, D. A 3D SPH–FE coupling for FSI problems and its application to tire hydroplaning simulations on rough ground. Comput. Methods Appl. Mech. Eng.; 2019; 355, pp. 558-590. [DOI: https://dx.doi.org/10.1016/j.cma.2019.06.033]

49. Martínez-Estévez, I.; Tagliafierro, B.; El Rahi, J.; Domínguez, J.M.; Crespo, A.J.C.; Troch, P.; Gómez-Gesteira, M. Coupling an SPH-based solver with an FEA structural solver to simulate free surface flows interacting with flexible structures. Comput. Methods Appl. Mech. Eng.; 2023; 410, 115989. [DOI: https://dx.doi.org/10.1016/j.cma.2023.115989]

50. El Rahi, J.; Martínez-Estévez, I.; Tagliafierro, B.; Domínguez, J.M.; Crespo, A.J.C.; Stratigaki, V.; Suzuki, T.; Troch, P. Numerical investigation of wave-induced flexible vegetation dynamics in 3D using a coupling between DualSPHysics and the FEA module of Project Chrono. Ocean Eng.; 2023; 285, 115227. [DOI: https://dx.doi.org/10.1016/j.oceaneng.2023.115227]

51. El Rahi, J.; Martínez-Estévez, I.; Almeida Reis, R.; Tagliafierro, B.; Domínguez, J.M.; Crespo, A.J.C.; Stratigaki, V.; Suzuki, T.; Troch, P. Exploring wave–vegetation interaction at stem scale: Analysis of the coupled flow–structure interactions using the SPH-based DualSPHysics code and the FEA module of Chrono. J. Mar. Sci. Eng.; 2024; 12, 1120. [DOI: https://dx.doi.org/10.3390/jmse12071120]

52. Domínguez, J.M.; Crespo, A.J.C.; Hall, M.; Altomare, C.; Wu, M.; Stratigaki, V.; Troch, P.; Cappietti, L.; Gómez-Gesteira, M. SPH simulation of floating structures with moorings. Coast. Eng.; 2019; 153, 103560. [DOI: https://dx.doi.org/10.1016/j.coastaleng.2019.103560]

53. Domínguez, J.M.; Fourtakas, G.; Altomare, C.; Canelas, R.B.; Tafuni, A.; García-Feal, O.; Martínez-Estévez, I.; Mokos, A.; Vacondio, R.; Crespo, A.J.C. . DualSPHysics: From fluid dynamics to multiphysics problems. Comput. Part. Mech.; 2022; 9, pp. 867-895. [DOI: https://dx.doi.org/10.1007/s40571-021-00404-2]

54. Tasora, A.; Serban, R.; Mazhar, H.; Pazouki, A.; Melanz, D.; Fleischmann, J.; Taylor, M.; Sugiyama, H.; Negrut, D. Chrono: An open source multi–physics dynamics engine. High Performance Computing in Science and Engineering; Springer International Publishing: Soláň, Czech Republic, 2016; pp. 19-49. [DOI: https://dx.doi.org/10.1007/978-3-319-40361-8_2]

55. Martínez-Estévez, I.; Domínguez, J.M.; Tagliafierro, B.; Canelas, R.B.; García-Feal, O.; Crespo, A.J.C.; Gómez-Gesteira, M. Coupling of an SPH-based solver with a multiphysics library. Comput. Phys. Commun.; 2023; 283, 108581. [DOI: https://dx.doi.org/10.1016/j.cpc.2022.108581]

56. Martínez-Estévez, I. MoorDyn+. Available online: https://github.com/imestevez/MoorDynPlus (accessed on 1 October 2025).

57. Hall, M.; Goupee, A. Validation of a lumped-mass mooring line model with DeepCwind semisubmersible model test data. Ocean Eng.; 2015; 104, pp. 590-603. [DOI: https://dx.doi.org/10.1016/j.oceaneng.2015.05.035]

58. Wendland, H. Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree. Adv. Comput. Math.; 1995; 4, pp. 389-396. [DOI: https://dx.doi.org/10.1007/BF02123482]

59. Batchelor, G.K. An Introduction to Fluid Dynamics; Cambridge University Press: Cambridge, UK, 2000; [DOI: https://dx.doi.org/10.1017/CBO9780511800955]

60. Fourtakas, G.; Dominguez, J.M.; Vacondio, R.; Rogers, B.D. Local Uniform Stencil (LUST) boundary condition for arbitrary 3-D boundaries in parallel smoothed particle hydrodynamics (SPH) models. Comput. Fluids; 2019; 190, pp. 346-361. [DOI: https://dx.doi.org/10.1016/j.compfluid.2019.06.009]

61. Oger, G.; Doring, M.; Alessandrini, B.; Ferrant, P. An improved SPH method: Towards higher order convergence. J. Comput. Phys.; 2007; 225, pp. 1472-1492. [DOI: https://dx.doi.org/10.1016/j.jcp.2007.01.039]

62. Michel, J.; Colagrossi, A.; Antuono, M.; Marrone, S. A regularized high-order diffusive smoothed particle hydrodynamics scheme without tensile instability. Phys. Fluids; 2023; 35, 103604. [DOI: https://dx.doi.org/10.1063/5.0165036]

63. Bonet, J.; Lok, T.-S.L. Variational and momentum preservation aspects of smooth particle hydrodynamic formulations. Comput. Methods Appl. Mech. Eng.; 1999; 180, pp. 97-115. [DOI: https://dx.doi.org/10.1016/S0045-7825(99)00051-1]

64. Marrone, S.; Colagrossi, A.; Le Touzé, D.; Graziani, G. Fast free-surface detection and level-set function definition in SPH solvers. J. Comput. Phys.; 2010; 229, pp. 3652-3663. [DOI: https://dx.doi.org/10.1016/j.jcp.2010.01.019]

65. Khayyer, A.; Gotoh, H.; Shimizu, Y. Comparative study on accuracy and conservation properties of two particle regularization schemes and proposal of an optimized particle shifting scheme in ISPH context. J. Comput. Phys.; 2017; 332, pp. 236-256. [DOI: https://dx.doi.org/10.1016/j.jcp.2016.12.005]

66. Altomare, C.; Crespo, A.J.C.; Domínguez, J.M.; Gómez-Gesteira, M.; Suzuki, T.; Verwaest, T. Applicability of Smoothed Particle Hydrodynamics for estimation of sea wave impact on coastal structures. Coast. Eng.; 2015; 96, pp. 1-12. [DOI: https://dx.doi.org/10.1016/j.coastaleng.2014.11.001]

67. Oger, G.; Marrone, S.; Le Touzé, D.; de Leffe, M. SPH accuracy improvement through the combination of a quasi-Lagrangian shifting transport velocity and consistent ALE formalisms. J. Comput. Phys.; 2016; 313, pp. 76-98. [DOI: https://dx.doi.org/10.1016/j.jcp.2016.02.039]

68. Lind, S.J.; Xu, R.; Stansby, P.K.; Rogers, B.D. Incompressible smoothed particle hydrodynamics for free-surface flows: A generalised diffusion-based algorithm for stability and validations for impulsive flows and propagating waves. J. Comput. Phys.; 2012; 231, pp. 1499-1523. [DOI: https://dx.doi.org/10.1016/j.jcp.2011.10.027]

69. English, A.; Domínguez, J.M.; Vacondio, R.; Crespo, A.J.C.; Stansby, P.K.; Lind, S.J.; Chiapponi, L.; Gómez-Gesteira, M. Modified dynamic boundary conditions (mDBC) for general-purpose smoothed particle hydrodynamics (SPH): Application to tank sloshing, dam break and fish pass problems. Comput. Part. Mech.; 2022; 9, pp. 1-15. [DOI: https://dx.doi.org/10.1007/s40571-021-00403-3]

70. Aalborg University. AwaSys. Available online: https://www.hydrosoft.civil.aau.dk/awasys/ (accessed on 1 October 2025).

71. Altomare, C.; Domínguez, J.M.; Crespo, A.J.C.; González-Cao, J.; Suzuki, T.; Gómez-Gesteira, M.; Troch, P. Long-crested wave generation and absorption for SPH-based DualSPHysics model. Coast. Eng.; 2017; 127, pp. 37-54. [DOI: https://dx.doi.org/10.1016/j.coastaleng.2017.06.004]

72. Leimkuhler, B.; Matthews, C. Molecular Dynamics: With Deterministic and Stochastic Numerical Methods; Springer: Cham, Switzerland, 2015; [DOI: https://dx.doi.org/10.1007/978-3-319-16375-8]

73. Monaghan, J.J.; Kos, A. Solitary waves on a Cretan beach. J. Waterw. Port Coast. Ocean Eng.; 1999; 125, pp. 145-155. [DOI: https://dx.doi.org/10.1061/(ASCE)0733-950X(1999)125:3(145)]

74. Tafuni, A.; Domínguez, J.M.; Vacondio, R.; Crespo, A.J.C. A versatile algorithm for the treatment of open boundary conditions in smoothed particle hydrodynamics GPU models. Comput. Methods Appl. Mech. Eng.; 2018; 342, pp. 604-624. [DOI: https://dx.doi.org/10.1016/j.cma.2018.08.004]

75. Liu, M.B.; Liu, G.R. Restoring particle consistency in smoothed particle hydrodynamics. Appl. Numer. Math.; 2006; 56, pp. 19-36. [DOI: https://dx.doi.org/10.1016/j.apnum.2005.02.012]

76. Meguro, K.; Tagel-Din, H. Applied element method for structural analysis: Theory and application for linear materials. Struct. Eng./Earthq. Eng. JSCE; 2000; 17, pp. 21-35. [DOI: https://dx.doi.org/10.2208/jscej.2000.647_31]

77. Meguro, K.; Tagel-Din, H. A new efficient technique for fracture analysis of structures. Bull. Earthq. Resist. Struct. Res. Cent.; 1997; 30, pp. 103-116.

78. Meguro, K.; Tagel-Din, H. Applied element method used for large displacement structural analysis. J. Nat. Disaster Sci.; 2002; 24, pp. 25-34.

79. Le Méhauté, B. An Introduction to Hydrodynamics and Water Waves; Springer Science & Business Media: Berlin, Germany, 1976.

80. Schreier, S.; Jacobi, G. Measuring hydroelastic deformation of very flexible floating structures. WCFS2020; Piątek, Ł.; Lim, S.H.; Wang, C.M.; de Graaf-van Dinther, R. Springer Singapore: Singapore, 2022; pp. 347-371. [DOI: https://dx.doi.org/10.1007/978-981-16-2256-4_21]

81. Dyneema. Available online: https://www.dyneema.com/ (accessed on 1 October 2025).

82. Zago, V.; Schulze, L.J.; Bilotta, G.; Almashan, N.; Dalrymple, R.A. Overcoming excessive numerical dissipation in SPH modeling of water waves. Coast. Eng.; 2021; 170, 104018. [DOI: https://dx.doi.org/10.1016/j.coastaleng.2021.104018]

© 2025 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/). Notwithstanding the ProQuest Terms and Conditions, you may use this content in accordance with the terms of the License.

A Versatile SPH Approach for Modelling Very Flexible and Modularized Floating Structures in Moored Configurations

Content area

Abstract

Full text

1. Introduction

2. Numerical Methods

2.1. SPH Fluid Solver

2.1.1. Governing Equations

2.1.2. Variable Resolution

2.2. MBD Structural Solver

2.2.1. Governing Equations

2.2.2. SPH-MBD Coupling

2.3. Lumped-Mass Mooring Solver

2.3.1. Governing Equations

2.3.2. SPH-MoorDyn+ Coupling

2.4. Applied Element Method

3. Validation

3.1. Moored Array of Interconnected Floating Rigid Pontoons

3.1.1. Experimental Setup

3.1.2. Numerical Setup

3.1.3. Results

3.2. Moored Very Flexible Floating Structure

3.2.1. Experimental Setup

3.2.2. Numerical Setup

3.2.3. Results

4. Discussion

5. Conclusions