Application of the immersed-body method to simulate wave-structure interactions

This study aims at demonstrating the capability of the immersed-body method to simulate wave–structureinteractionsusinganon-linearfinite-elementmodel.Inthisapproach,theNavier–Stokes equations are solved on an extended mesh covering the whole computational domain (i.e. fluids and structure). The structure is identified on the extended mesh through a nonzero solid-concentration field, which is obtained by conservatively mapping the mesh discretising the structure onto the extended mesh. A penalty term relaxes the fluid and structural velocities to one another in the regions covered by the structure. The paper is novel in that it combines the immersed-body method with wave modelling and mesh adaptivity. The focus of the paper is therefore on demonstrating the capability of this new methodology in reproducing well-established test cases, rather than investigating new physical phenomena in wave–structure interactions. Two cases are considered for a bottom-mounted pile. First, thepileisplacedinanumericalwavetank,wherepropagatingwavesaremodelledthroughafree-surface boundary condition. For regular and irregular waves, it is shown that the wave dynamics are accurately modelled by the computational fluid dynamics model and only small discrepancies are observed in the close vicinity of the structure. Second, the structure is subjected to a dam-break wave impact obtained by removing a barrier between air and water. In that case, an additional advection equation is solved for a fluid-concentration field that tracks the evolution of the air–water interface. It is shown that the load associated with the wave impact on the structure compares well with existing numerical and experimental data.


Introduction
The accurate computation of wave loads is important in offshore engineering, for example, to optimally design oil and gas platforms or coastal defence structures. In recent years, there has also been an increased interest in offshore renewableenergy devices, which can be either mounted on the sea bed (fixed devices) or moored to it (floating devices). In this context, an accurate prediction of the wave loading is vital to ensure that offshore renewable-energy devices are economically viable and can withstand rough sea conditions. Computational models can assist in the design of offshore renewable-energy devices by analysing several different configurations, while limiting expensive laboratory or onsite testing. The hydrodynamic behaviour of such devices is however complex due to: (i) the interactions between extreme waves and structures (either fixed or floating), and (ii) the mutual interactions between fluids and structures if the latter floats. The method presented in this paper has been developed to tackle both aspects. Only fixed structures are considered in this manuscript.
Various methods exist to numerically model fluids interacting with a structure. One possibility consists in solving the Navier-Stokes equations on a mesh excluding the structure (defined-body method). This is typically done if the structure is fixed. When the structure moves in the fluid domain, re-meshing is necessary. This is computationally expensive and might also yield highly-distorted grids. Another approach, which does not require re-meshing, consists in meshing the entire domain (containing both fluids and structure) and modelling the effect of the structure through surface or body forces. This methodology underpins the so-called immersed-boundary method, proposed by Peskin [1] in the context of cardiac mechanics. Reviews of immersed-boundary type methods are available in the literature [2][3][4]. This method has been successfully applied to many flows of technological relevance, partly because it is very efficient in dealing with the presence of complex boundaries [5,6].
The use of nonlinear computational fluid dynamics models in the context of wave-structure interactions is the exception rather than the rule. Various alternatives exist, such as empirical methods, physical testing, and simplified numerical approaches, as detailed by Folley et al. [7] in the context of wave energy converter arrays. If potential flow is assumed, the scattering of plane waves induced by a bottom-mounted cylinder, in the linear diffraction regime, is described by a Laplace equation for the velocity potential and a set of boundary conditions. The problem has an analytical solution [8]. It is also commonly solved numerically using the so-called boundary element method, which works on either a linearised formulation of the problem (linear potential models) or the non-linear formulation (non-linear potential models). Linear approaches operate in the frequency domain and are therefore not computationally demanding. Both their low computational demand and limited complexity make them ideally suited to industry-standard applications. In contrast, non-linear approaches operate in the time domain, so that both transient phenomena and non-linear wave-structure interactions can be addressed. Therefore, non-linear potential models are increasingly used for computing extreme loads on fixed and moving structures. Computational fluid dynamics models are attractive when potential flow assumptions cannot be made, because they can account for viscous and rotational effects through the resolution of the Navier-Stokes equations. This is particularly needed when the structure is subjected to extreme waves, for which viscous effects and air entrainment cannot be neglected. However, the accurate prediction of wave propagation using computational fluid dynamics models is difficult, even in the absence of a structure. This is due to the inherent energy dissipation introduced by the discretisation schemes [9]. The ability of the finite-element model 'Fluidity' [10,11] to simulate linear and moderately non-linear gravity waves was demonstrated in two-and three-dimensional numerical wave tanks [12,13]. Good agreement was found between the simulation results and the theoretical predictions in the inviscid regime.
The purpose of this work is to analyse the capability of 'Fluidity' in the simulation of waves interacting with a structure. The structure is modelled as a volumetric body force immersed in the fluid domain [14]. The so-called immersed-body approach is a continuous forcing approach based on a penalty term. It relies on two distinct meshes: one covering the entire domain (extended mesh), and the other discretising solely the structure (solid mesh). The region occupied by the structure is identified through a solidconcentration field, which is computed by projecting a unitary field from solid to extended mesh. One of the advantages of this dualmesh approach is that each mesh can be used by a different finiteelement model when both the fluid and the structural dynamics are of interest, as shown in [14]. The novelty of this paper stems from the combination of the immersed-body approach with wave modelling and adaptive re-meshing. Since this methodology has never been applied to wave-structure interactions, the focus of the paper is on well-established test cases rather than on studying new physical phenomena. The paper is organised as follows. The general principle of the immersed-body method and the governing equations are derived in Section 2. Section 3 details the methods used in this study to model waves and Section 4 focuses on the procedure used to compute the load on the structure. Section 5 explains how the fluid-dynamics model can modify the extended mesh dynamically in time, in order to refine the resolution around certain flow features. Dynamic mesh adaptivity is particularly attractive to reduce the computational cost of fluid dynamics simulations. Finally, results are shown in Section 6 for three cases. First, a cylindrical pile is subjected to a regular train of smallamplitude gravity waves in the inviscid regime. In this context, linear diffraction theory [8] is used as reference solution. Second, the same cylindrical pile is subjected to an irregular focused wave group. In considering this case, the surface elevation in the absence of the structure is first compared to the secondorder irregular wave solution by [15]. Third, a square cylinder is subjected to a dam-break water wave in the viscous regime. Results are compared with numerical and experimental results from the literature [16]. Conclusions are drawn in Section 7.

Underlying equations for modelling fluid-structure interactions
The fluid/marine dynamics model 'Fluidity' is a finite-element open-source (LPGL) numerical tool which solves the non-hydrostatic Navier-Stokes equations on unstructured meshes [17,10,18]. In a purely hydrodynamics problem, the equations of motion of the fluids are given by where u f is the fluid velocity, p is the pressure field, ρ is the fluid density, µ is the dynamic viscosity of the fluid, S is the deviatoric part of the stress tensor, and B represents other source terms (such as the gravitational force). When a structure is immersed in the fluid domain, the fluid-structure interactions can be modelled in two different ways: either by excluding the structure from the computational mesh and solving only for the fluid dynamics in the regions covered by the fluid (defined-body method), or by solving the equations of motion in an extended domain that covers both fluid and structure (immersed-body method [14]). The definedbody method is typically used for fixed structures. If the structure moves, however, it becomes computationally expensive because re-meshing is necessary to track the structure's motion. It can also lead to highly-distorted grids. In this context, the immersedbody method is attractive because re-meshing is not necessary. Instead, the regions V s covered by the structure are filled with the surrounding fluid V f so that one computational mesh spreads across the whole domain (i.e. V = V f ∪ V s ). At the location of the structure, a penalty term is used to relax fluid and structural velocities to one another. The momentum equations for the fluids are therefore modified in two ways. First, they have to be solved for a monolithic velocity, which is defined over the extended domain V , i.e.
where u s is the solid velocity, and the concentration fields are defined as α f = V f /V for the fluid and α s = V s /V for the solid.
By definition, α f + α s = 1 is satisfied in the whole domain.
Second, a penalty force F is added to relax the monolithic and solid velocities to one another. Therefore, the momentum balance on the extended mesh becomes where the penalty force F weakly ensures that the monolithic and solid velocities equal one another in the structure. The penalty force is thus expressed as where the relaxation factor β = ρ f / t dictates how fast the monolithic and solid velocities equal one another. In this study, it is assumed that β is driven by inertial effects rather than viscous effects. This assumes that ρ/ t ≫ µ/l 2 e , where t is the computational time step and l e is the local edge length in the mesh. The penalty term is nonzero only in the regions covered by the solids. These regions are identified by a nonzero solid-concentration field α s on the extended mesh. The solidconcentration field is computed by conservatively mapping the solid mesh onto the extended mesh [14]. The mapping is achieved by projecting a unitary field from solid to extended mesh using a Galerkin projection, so that is satisfied at a discrete level (F denoting the extended mesh and S being the solid mesh). Importantly, the monolithic velocity directly results from the resolution of Eqs. (4)-(5) over the extended computational mesh, rather than being derived from Eq. (3). The equations are discretised spatially using finite elements and temporally using a Crank-Nicolson scheme, as already reported in [14,11].

Wave modelling
There are two ways of modelling waves in 'Fluidity'. The first approach consists in defining a computational domain that contains only water and representing waves as a free-surface boundary condition [12]. In this framework, Eq. (5) is solved under the Boussinesq approximation. Therefore, the density at a position x can be written as ρ( As a result, the total density can be replaced by ρ 0 in all the terms of Eq. (5), except in the gravitational force. The Navier-Stokes equations then become ∇ · u = 0, whereû accounts for the velocity of a mesh moving with the free-surface, η is the free-surface elevation above a reference level z = 0, p ′ is the perturbation pressure, and g is the gravitational acceleration pointing in the direction k. The velocity field at the free-surface is related to the free-surface elevation through the kinematic boundary condition, which imposes that the velocity of a water particle normal to the surface is equal to the speed of the surface in that direction [19]. The algorithm used for computing the free-surface elevation is detailed in [20]. This method is used in Sections 6.1 and 6.2, where waves propagate in a numerical wave tank. In this context, a mixed finite-element discretisation method is used, where discontinuous linear polynomials are used for the velocity field (i.e. a P 1DG discretisation) and continuous piecewise quadratic polynomials are used for the pressure field (i.e. a P 2 discretisation) [21,22]. The second approach consists in defining a computational domain that contains both air and water. Thus, the governing equations are solved for each fluid, assuming that the fluid phases are immiscible. The Navier-Stokes equations for the fluid element i are expressed as (11) where ζ i denotes the fluid-concentration field and varies between 0 and 1. An additional advection equation is solved for ζ i , in order to track the evolution of the air-water interface, i.e.
The interface-tracking algorithm of 'Fluidity' is extensively described by Wilson [23]. This approach is used in Section 6.3, where the incoming wave is generated through the collapsing of a water column. In that case, the advection equation for the fluidconcentration field is solved on a dual control-volume mesh, for a piecewise-constant representation of the fluid-concentration field. Consistent discretisations for the momentum and material advection steps further ensure conservation of ζ i . Additionally, a high-order flux reconstruction in combination with a flux limiting scheme is used to accurately track the air-water interface on unstructured meshes. This provides a simple alternative to schemes that require an explicit interface reconstruction, and particle based methods [23]. The advective fluxes are limited using the HyperC approach, and the Bassi-Rebay discretisation scheme is used for the diffusion term [24,18]. The velocity field is piecewise constant over the elements (i.e. a P 0 discretisation), while the pressure field is piecewise constant over control volumes (CV) and its degrees of freedom are stored in exactly the same locations as P 1 (i.e. a P 1CV discretisation). The advantage of the P 0 − P 1CV discretisation pair is that the advective velocity will be exactly divergence-free in the advection equation for P 1CV tracers, as the continuity equation is tested with P 1CV test functions.

Load computation
As mentioned in the introduction, an accurate prediction of the wave loading is crucial in offshore engineering applications. In order to compute the load on the structure, the pressure and velocity fields are first projected from the extended mesh to the solid mesh using a Galerkin projection. Given a donor mesh D and a target mesh T on a domain Ω, a Galerkin projection on a field q from D to T ensures that by minimising the L 2 norm of the interpolation error. These projections are explained in detail in other publications by the authors [14,25], where the forces are transferred from a fluiddynamics solver to a solid-dynamics solver (and vice-versa) in the context of fluid-structure interactions with moving (rigid and deformable) structures. In a discrete form, Eq. (13) involves solving a linear system of equations containing a mass matrix based on the basis functions of the target mesh and a mixed mass matrix formed of the basis functions of the target and donor meshes. In this study, the intersections between both meshes are further identified using an advancing front algorithm [26].
Once the projections are completed, the pressure gradient and divergence of the Reynolds stresses are computed on the solid mesh. The total load on the structure is obtained by integrating their sum over the solid mesh, i.e.
According to Gauss' theorem, this is equivalent to integrating pressure and Reynolds stresses over the surface of the structure. Following this procedure, the accuracy of the load calculation depends on the resolution of the solid mesh. A detailed analysis on the mesh requirements is beyond the scope of this study. Here, the resolution of the solid mesh is chosen such that it yields an accurate determination of the load.

Mesh adaptivity
'Fluidity' has the capability of optimising the mesh in time, in order to refine the spatial resolution around certain flow features. For example, if vortices develop in the wake of a structure, the mesh can be periodically re-generated to focus the resolution around the shed vortices and coarsen it elsewhere. Details on the mesh-adaptivity procedure can be found in [27]. In this work, a hr-adaptive method is used where both the mesh vertices and their connectivity are changed throughout the simulation. In this context, the desired geometric properties of the new mesh are described by a metric, which is a symmetric positive-definite tensor field. The metric is computed from the Hessian H of the solution field and a user-defined weight ϵ, i.e.
where n is the dimension of the space and the modified Hessian In Eq. (16) According to Eq. (15), the mesh is refined for high curvatures in the chosen field (i.e. large eigenvalues) or small values of the field weight ϵ. Conversely, the mesh is coarsened in regions of low curvature in the field or if the field weight is increased. In this study, mesh adaptivity is controlled solely by the target value ϵ of the L 2 -norm of the solution's interpolation error. Mesh adaptivity is used in Section 6.3.

Small-amplitude waves interacting with a cylindrical pile
The ability of 'Fluidity' to simulate small-amplitude regular waves propagating in a numerical wave tank (without structure) has already been demonstrated [12,13]. In this section, a cylindrical pile is mounted in the numerical wave tank to investigate wave-structure interactions. Fig. 1 shows a sketch of the domain. The centre of the cylindrical pile is placed at a distance of 42D from the tank outlet, and 17D from the inlet and sides (D being the pile diameter). The tank depth is h = 2.4D. The horizontal, lateral and vertical directions are denoted by x, y, and z respectively. Regular small-amplitude waves are generated by setting the horizontal velocity at the inlet (x = 0) to the linear wave solution, where a is the wave amplitude and the wavenumber k is given by the associated dispersion equation, In Eq. (19), the wave frequency is ω = 2π /T , where T is the wave period. In this paper, the simulations are performed at full scale, with a pile length that is typical for fixed or floating wind turbines.
Therefore, the wave period is chosen to be equal to T = 10 s, as typically used for the design of such structures. Small-amplitude waves are considered for two different values of the wave steepness ak: a very small amplitude wave with ak = 0.001 and a finite, yet small, amplitude wave with ak = 0.01.
In the latter, the linear diffraction theory is still valid even if the deformations of the free-surface become significant compared to the structure diameter. Defining the deep-water wavelength as λ 0 = 2π g/ω 2 , the present non-dimensional water depth in all cases is equal to h/λ 0 = 0.45, which is considered as an intermediate water depth [12]. At the outlet of the domain (x = 59D), the velocity components are left free and the non-hydrostatic part of the pressure field is naturally set to zero. An absorption layer is also used in the region 51D ≤ x ≤ 59D to avoid spurious wave reflections. To this end, the artificial absorption term σ u is added to the right-hand-side of the Navier-Stokes equations, where withx = (x − L 0 )/L, and L = 8D is the length of the absorption layer that commences at x = L 0 = 51D [28]. A no-normal flow condition is used at the bottom and lateral sides of the tank, while a combined pressure and free-surface kinematic boundary condition is prescribed for the top surface [12]. An unstructured mesh is generated in the x-y plane and extruded vertically into tetrahedra using seven uniformly-spaced layers. The typical edge length of the mesh elements in the horizontal plane varies from l e = D/2 at the external boundaries to l e = D/10 around the cylinder. Fig. 2 shows the extended The flow is considered inviscid to enable direct validation against linear diffraction theory [8]. The present numerical results are shown for modelling the pile using two different approaches: either the computational mesh excludes the region occupied by the pile (defined-body method) or the mesh covers the entire wave tank and a separate solid mesh discretises the three-dimensional pile (immersed-body method). In the latter case, the position of the pile on the extended mesh is represented through a solidconcentration field, as described in Section 2. Figs. 3 and 4 show the evolution of the free-surface elevation in the horizontal direction, at the pile centreline and t = 15.5T , and for different values of wave steepness. Symbols represent the numerical result using the defined-body method (left) and the immersed-body method (right), while the continuous line shows the analytical solution by [8]. The grey area further illustrates the location of the pile. The overall agreement between the computational fluid dynamics results and the linear diffraction theory is very good. However, local discrepancies are observed in the vicinity of the pile. In particular, the value of the free-surface elevation in the definedbody case at the front and rear of the pile is within a maximum of 3% of the theoretical solution for both values of wave steepness investigated. In the immersed-body case, the discrepancies vary between a 7% underestimation of the value at the front and a 15% overshoot at the rear. These percentages decrease to 3% and 11% at t = T /100. Similar observations can be made on the horizontal velocity field u x , whose evolution in the horizontal direction is shown in Figs. 5 and 6 for the two values of the wave steepness. A snapshot of the horizontal velocity profile is also shown in Fig. 7 for ak = 0.001 and a depth of z = −h/7. A possible explanation for larger discrepancies obtained with the immersedbody approach in comparison to the defined-body approach is that the representation of the immersed pile by the solid-concentration field is not sharp. The solid-concentration field linearly varies between 0 (at the physical location of the pile boundary) and 1 (at the first mesh node inside the pile). Although the pile is reasonably well represented on the computational mesh, the value of the freesurface at the first mesh node inside the pile is nonzero. These nonphysical values could explain the presence of slight overshoots and undershoots of the free-surface elevation around the pile, when it is modelled as an immersed body. Another explanation is that the penalty forcing used in the immersed-body approach tends to impose a no-slip boundary condition as physically encountered on the surface of the structure. Discrepancies at the fluid-structure interface could therefore be observed when the computation is inviscid. Finally, the fact that there is no major differences between the results obtained for the very small amplitude wave (Fig. 3) and the finite amplitude wave (Fig. 4) is encouraging. It confirms that wave energy is conserved in the computational fluid dynamics model, as already observed in the absence of a structure [12].

Nonlinear waves interacting with a cylindrical pile
The regular waves considered in the previous section are not very realistic because, in practise, sea states must be considered as irregular. In this section, the computational domain is similar to the numerical wave tank considered in the previous section, except that the wave maker generates relatively steep transient wave groups; these transient wave groups being achieved through the focusing of multiple frequency components. Also no absorption layer is used because the analysis concerns the results close to the time at which the focused event takes place, i.e. when the waves are not reflected by the domain boundaries. The focused wave events are chosen such that the maximum wave amplitude occurs at a distance x f = 10 h from the input boundary. Three values are considered for the sum of focused event amplitude: A sum k p = 0.018, A sum k p = 0.09, and A sum k p = 0.18, where k p is the wavenumber corresponding to the peak in the JONSWAP spectrum [29] and A sum corresponds to the linear amplitude sum of all wave components. The latter has a spectral density function given by where  Additionally, ω p = 2π /T p is the circular peak frequency, T p is the peak period, α is the energy scale parameter, and γ is the peak enhancement factor. For the purpose of the present work, these parameters were chosen as γ = 2.5 and T p = 16 s, with α being adjusted to achieve a particular amplitude sum A sum . The wavemaker input kinematics were computed using the secondorder solution by [15]. Since it is numerically challenging to have a stable solution for large values of the wave steepness when the flow is inviscid, a uniform kinematic viscosity equal to ν = 10 −3 m 2 s −1 is used in the present simulations as viscous damping. Ongoing work is investigating different stabilisation methods that could minimise this viscous dissipation.
First, the numerical results are shown in the absence of a pile. In that case, the tank width is reduced to a value of w = 0.34D, because the waves are uni-directional and undisturbed by diffraction effects. The tank dimensions in the streamwise and depth-normal directions are 51D and 2.4D, respectively (D being the pile diameter as defined in the previous section). The edge Good overall agreement is obtained between the numerical results and the second-order solution at A sum k p = 0.018 (Fig. 8) and A sum k p = 0.09 (Fig. 9). The slight attenuation of the maximum crest elevation is believed to be due to the viscous damping term. At larger amplitudes of the focused event (Fig. 10), the computational fluid dynamics results exhibit a downstream shift of the focused location (compared to the second-order solution) and an asymmetry in the minima of the free-surface elevation.
This is expected because, at A sum k p = 0.18, harmonics higher than the second order become non negligible as demonstrated by [30]. Therefore, these results show that the computational fluid dynamics model is capable of capturing high-order harmonics in the free-surface elevation. Perhaps more importantly, the ability of modelling this freesurface effect can now be combined with the immersed-body method. To achieve this, the same pile as used in Section 6.1 is mounted in the tank at a distance of 10.   percentages further increase to reach 28% and 43%; the difference between the two cases clearly emphasising the importance of nonlinearity. Finally, Fig. 12 highlights the spatial evolution of the non-  dimensionalised surface elevation for A sum k p = 0.09 and t = 7.5T p using the defined-body method (left) and immersed-body method (right). Note that the nonzero value of the free-surface elevation inside the pile in the defined-body case is an artefact of the colour plot. The same diagnostics are shown at t = 8.75T p by Fig. 13.
These figures can be further compared to Fig. 14, which shows the same diagnostics in the absence of pile. The diffraction effects due to the presence of the pile are particularly apparent at t = 8.75T p .

Dam-break wave impact on a square cylinder
This test case differs from the previous ones in that the dynamics of two fluids (air and water) are considered. The initial domain configuration is shown by Fig. 15 (left). A square pile of width D and height 6D is placed in a domain of size 13D × 5D × 6D.
The pile is centred at a distance 8D from the left boundary of the domain. The latter is delimited by solid boundaries, on which the normal component of the velocity field is set to zero. A water reservoir of height h 0 = 2.5D, length 3.5D and width 5D is and Re 2 = ρ 2 u p D/µ 2 = 9.6 · 10 6 , where u p is the peak value of the vertical velocity when the fluid impacts upon the pile. The air-water interface is tracked by solving an advection-diffusion equation for a fluid-concentration field, as explained in Section 3. Additionally, mesh adaptivity is used on the fluid-and solid-concentration fields with weights equal to ϵ α f = 0.09 and ϵ α s = 0.001. The minimum value of the mesh edge length is further fixed at D/50 in the adaptive algorithm. As a result, the adapted mesh has approximately 1.2 million nodes. Finally, the time step size is such that the Courant number is fixed at 0.1. In this work, it is important to highlight that the simulations are performed without subgrid scale model despite the fact that the Reynolds number is large. This means that the smallest turbulent structures are unresolved. The study of wave-structure interactions with turbulence modelling is the subject of ongoing research. Fig. 16 shows the time evolution of the force acting on the cylinder. Time and force are non-dimensionalised by √ h 0 /g and F ref = ρ 2 gh 2 0 D/2, respectively. The dashed and dashdotted lines highlight the results obtained with the defined-and immersed-body approach, respectively. The continuous line is the numerical result available in the literature [16], which considers the same setup. Finally, the dots show the envelope delimiting the experimental results from [16]. The present results agree well with the reference solutions, especially for predicting the force due to the incoming wave. At t ⋆ ≈ 7, the load on the cylinder drops and changes sign. This is because the incoming wave has hit the outer boundary of the domain yielding a reflected wave. The time at which the latter impacts upon the pile is slightly under-predicted in the present work. This could be due to two factors. First, the finite-element pair used is of relatively low order (a P 0 -P 1CV discretisation). As a result, the accuracy of the wave tracking tends to decrease as time evolves due to numerical dissipation. Second, the spatial resolution used at the outer boundary is limited because most of the mesh refinement takes place around the pile, as illustrated by Fig. 15 (right). This affects the accuracy with which the reflected wave is modelled. The detailed analysis of the flow features close to solid boundaries is the subject of ongoing work.

Conclusions
This paper focuses on the application of the immersed-body method to the simulation of wave-structure interactions. The immersed-body method is a versatile way to model fluid-structure interactions. It consists of immersing the structure in an extended fluid domain and weakly applying the boundary condition at the fluid-structure interface through a penalty force. In this paper, results obtained with the immersed-body method were systematically compared to those obtained by excluding the structure from the computational domain (defined-body approach). Two different numerical setups were considered. First, gravity waves were propagated in a numerical wave tank filled only with water (one-fluid problem). For regular waves, the paper showed that the overall wave diffraction behaviour was well reproduced by both the defined-and immersed-body approaches. The weakness of the immersed-body approach lay in the values of the wave elevation in the close vicinity of the structure, where overshoots and undershoots of the wave elevation are observed. These local effects are potentially of importance for the load computation, and therefore, require further analysis. However, they do not seem to significantly affect the wave dynamics away from the structure. Good agreement between defined-and immersed-body approaches was also obtained when irregular focused wave groups interacted with the pile. Second, a dam-break wave was generated by the collapse of a water column. In this context, an additional advection equation was solved for the air-water interface. The load associated with the wave impact on the structure was accurately modelled using both defined-and immersed-body approaches. A challenge in tracking fluid-fluid interfaces is however to reduce the numerical dissipation associated with low-order discretisation methods. Future work will focus on the in-depth analysis of flow features in the vicinity of the structure as well as incorporating high-order discretisation schemes or parametrisations.