# Simulations of a Partially Submerged Cylinder Subjected to Current and Waves

Wojciech Aniszewski<EMAIL_ADDRESS>Valentin Ageorges  Gaële Perret  Vincent
Moureau<EMAIL_ADDRESS>CORIA, CNRS UMR6614, Normandie Université, UNIROUEN,
INSA of Rouen, France  Laboratoire Ondes et Milieux Complexes (LOMC) – CNRS
UMR 6294 – Le Havre Normandy University, France

###### Abstract

This paper presents a new method allowing for a numerical study of the
interaction between a submerged solid structure and a current with surface
waves. The method involves a Navier-Stokes equation solver for a two-fluid
formulation with the Accurate Conservative Level-set (ACLS) method used for
interface tracking. The solver is generalized on unstructured grids, thus
enabling a boundary-fitted implementation of the FSI (Fluid-Solid
Interaction). A solid, vertical, partially submerged cylinder is examined as
it moves through the flume. The latter is partially filled with flowing water;
waves propagating on the liquid surface impact the cylinder. Simulation
results include examination of the free surface topology, air entrapment and
vorticity production in the cylinder wake, as well as the pressure
distribution on the solid. Numerical results are validated against an
experiment performed by the same authors.

###### keywords:

Wave generation, , fluid-structure interaction (FSI), simulation, Acurate
Conservative Level Set (ACLS)

††journal: International Journal of Multiphase Flow

##  1 Introduction

This paper presents a numerical study of the fluid-structure interaction in a
realistic geometry, validated by an accompanying experimental investigation.
The research presented here concerns a rigid structure in form of a partially
submerged, vertical cylinder, which interacts with a two-phase flow – namely a
current with surface waves. In the laboratory setup, this takes place inside a
$30$m-long flume. This type of flows occurs naturally in a number of coastal
engineering problems, such as oil rigs subjected to oceanic waves (Cleary and
Rudman, 2009), coastal installations subject to constant hydrodynamic loading,
ship hydrodynamics or construction of bridge pillars. Depending on the scale
of the problem, the studies in this broad field range from individual
structures causing structure blocking (Santo et al., 2018) to bridge damage by
hurricanes (Padgett et al., 2008).

Studies of fluid-structure interactions, while spanning a wide range of
spatio-temporal scales, dimensionless numbers and characteristic structure
geometries, share certain challenges. Examples of these challenges are found
in the realistic representation of both the liquid(s) and the solids in the
numerical simulation, or reliable and feasible measurement of the pressure
and/or velocity field in the experiment – especially for moving structures
(see e.g. Lobovský et al., 2014). Thus, regardless of the shape or character
of the structures or the nature of the flows, some of the problems will be
universal. For numerical simulations, these include e.g. the robustness of the
interface representation, especially near the solid. For example, the work of
Fondelli et al. (2015) presents an interesting study of the hydrodynamic
loading by means of both simulation and experiments. Namely, the authors are
focusing on the dam-break flow in which the bulk liquid, once freed from the
initial containment, impacts a solid, motionless block placed in its way on
the container bottom (and equipped with pressure sensors). The authors of
Fondelli et al. (2015) manage to reproduce the inertia-dominated stages of the
flow, while the subsequent sloshing stages prove more problematic, most likely
due to aforementioned shortcomings in the interface tracking scheme – which in
their case is the Volume-of-Fluid method with the adaptive refinement of a
structured mesh.

Ai et al. (2019) and Xie et al. (2020) present almost identical numerical
simulations of a 3D singular wave interacting with a vertical cylinder (in
both cases many more test cases are included for the presented codes). To
represent the interface, respectively the free-surface non-hydrostatic model
and the simplified algebraic version of the VOF scheme (see e.g. Aniszewski et
al., 2014) are used. The channel depth there is only $0.04$m, an order of
magnitude smaller than presented in this work. The authors of Ai et al. (2019)
impose the normal velocities of the incident periodic wave. Both Ai et al.
(2019) and Xie et al. (2020) include comparisons between simulations and
experiment for the solitary wave interacting with the cylinder. Interestingly
in the case of both papers, while the incident wave elevation and the
hydrodynamic load are predicted well, there are discrepancies in the cylinder
wake for both of the quantities. This is probably due to the deficiency of the
simple interface descriptions used. This suggests that it would not be trivial
to apply the simplistic interface treatment to the cases presented in our
study, which include liquid overturning and/or air bubbles entrained in the
cylinder wake. Other studies of submerged structures of varying shapes can be
found. For example, the authors of Arabi et al. (2019) study the interaction
of a solitary wave with a rectangular impervious block. Their findings –
including a strong correlation between the structure’s edges and the induced
vortical structures in the wake – concern also the mechanism of vorticity
production at the bottom of the partly submerged cylinder. A similar study to
the one presented here is found in the work of Corvaro et al. (2019), wherein
the authors use commercial software to study the wave-cylinder interactions at
lower Reynolds numbers.

Considering wave interactions with submerged structures, there exist a wide
array of analytical (Taylor and Wu, 1997; Davis, 1974) and numerical (e.g.
Brocklehurst et al., 2011; Han et al., 2015) works aiming at various
structures - such as bridges, bridge decks and pillars, including flexible
ones. For example, the work Qu et al. (2021) presents a numerical
investigation into the wave interactions with a coastal structure. While their
focus is on the ”extreme” waves (large amplitude oceanic waves originating
from spontaneous focusing), the validation cases presented therein fit in the
same flow regimes as presented in our work. Notably, the wave impact loads on
a rectangular 2D cylinder are studied with results comparing reasonably to the
experimental data. Moreover, the work of Qu et al. (2021) contains a numerical
wave generation technique at the channel inlet similar to this work (see
Sect.2.3). However, the work of Qu et al. (2021) – or comparable studies such
as Hu et al. (2016) – exhibits possible deficiencies in the interface
description or simply lack of spatial resolution. This in turn could have made
it harder for their authors to trace phenomena such as near-structure air
entrapment or drag reduction in the wake. This problem is tackled in our work
thanks to the unstructured, boundary-fitted meshes with variable resolution
and a massively parallel code.

The doctoral thesis one of these authors (Ageorges, 2019) and the companion
publications (Ageorges et al., 2019, 2021) contain a study of the air
entrapment in the flow around a partially submerged cylinder. The
investigations there are both experimental and numerical using the YALES2 code
(Moureau et al., 2011) used in this work as well, for a wide range of Reynolds
and Froude numbers. These works do not involve water waves: surface of the
liquid is at rest before impacting the (moving) solid cylinder. The work
identified two types of air entrapment (in a cavity along the cylinder wall
and in the cylinder wake), and examined the critical velocities at which this
phenomenon appears. Considering the cylinder drag, the work of Ageorges et al.
(2019) confirmed that the drag coefficient is generally smaller in multi- than
in single-phase flow, an effect made even stronger by the air entrapment. The
works in Ageorges et al. (2019, 2021) are comparable to certain studies
published using structured meshes (see e.g. Elhanafi et al. (2021)) but avoids
inherent numerical problems such as the solid modelling on rectangular grids,
or limited resolution of the CFD simulations. The work presented here aims to
continue the studies presented in Ageorges et al. (2019, 2021), by first
refining the numerical methodology, and subsequently including the modelling
of the incident waves in the fluid-structure interaction study. Unlike the
previous works, here we focus on a single flow configuration – corresponding
to an intermediate value of the Reynolds number, intermediate wave amplitudes
and low intensity of air entrapment in the cylinder wake – and present a
throughout comparison of the experimental and numerical approaches.

#### Structure of the paper

This paper is a part of a wider research project aimed at studying the liquid-
structure interactions in varying configurations of hydrodynamic loading in
single and multi-phase flows. Since our work involves both the experimental
facility as well as the simulation framework, it is necessary to restrict the
presentation here to its single aspect, for the sake of clarity. Thus, this
paper is focused on presenting the proposed numerical methodology and the
evidence of its validity. This dictates the following structure of the paper.
In the Section 2 we introduce the numerical solver YALES2, describing the
actual manner of solving the governing equations (2.1), generation of the
computational meshes (2.5) the wave generators (2.3) and the post-processing
tools (2.4). Subsequently, we present a concise description of the experiment
(Section 3 along with references to more publications consecrated to it). In
the Section 4 we focus our attention on the results of the simulations, which
we study in details including the macroscopic flow evolution (4.2) or the
spectral analysis of certain signals (4.5). The evolution of the pressure
force acting on the cylinder is described in Section 4.3 for both the
experiment and simulations. This is followed by a more detailed analysis of
the interfacial waves (e.g. the reflection coefficient in 4.6), analysis of
the flow statistics (4.7), after which the paper concludes.

##  2 Numerical methodology

###  2.1 Code Characteristics

YALES2111Yet Another Large-Eddy Simulation Solver is a massively parallel,
finite volume - based, structured and unstructured grid code, which began its
existence as a combustion simulation tool (Moureau et al., 2011) for
realistic, industrial geometries. It has since matured into a fully-fledged
multi-physics code including nearly $40$ solvers. The code has been created at
CORIA and is currently maintained by CORIA and a large global community of
contributors from academia and industry222More information can be found on the
YALES2 public website: https://www.coria-cfd.fr/index.php/YALES2.. The code is
written in Fortran, and employs – within the F90 and F2008 language standards
– data structures and methods characteristic for object-oriented programming.
Recent applications of YALES2 include the two-phase atomizing flow studies
(Janodet et al., 2022), new approaches to simulating wind turbines (Benard et
al., 2018), simulations of the Covid-19 viral transmissions (Dubief et al.,
2020) or proposing new numerical framework for high order approximation on
finite volumes (Bernard et al., 2020).

To characterize the code briefly, we note that YALES2 is distinguished mainly
by the internal data structures, created in the bottom-up fashion with massive
parallelism in mind. The main grid object can link to a structured or
unstructured mesh, with the latter being more optimized (and the former used
mostly for testing and benchmarks in academic geometries). Thus, in realistic
geometries, boundary fitted approach is used to model the domain boundaries.
When working in parallel, most widespread approach to the CPU load
distribution problem is Domain Decomposition (DD, see e.g. Chen et al., 2017),
i.e. splitting the domain spatially and redistributing the grid nodes between
the computing threads in a dynamical or static fashion. In YALES2, double-
domain decomposition (DDD) is used, meaning not only thread (CPU) domains
exist, but each of them has sub-domains called element groups. An additional
communication layer (internal communicator) exists in YALES2 on top of the
external one (MPI), and the former permits information exchange between
element groups and synchronization. This is justified on large distributed
systems such as super-computers because while the cost of MPI communication
remains similar in YALES2 to that using DD, its application of DDD greatly
diminishes execution slowdowns due to memory cache misses and memory access on
each of the computing cores. An additional benefit of DDD is that on each
computing core, certain operations become possible without any MPI
communications with neighbouring cores. An example of such an operation is
generation of a coarse mesh e.g. for conditioning in the solution of the
Poisson equation. To maintain optimal partitioning on all levels, the external
METIS library (Karypis, 2003) is used. YALES2 is able to import simulation
geometry in many formats (i.e. Ansys, HDF5, Gambit, GMesh) and is able to
optimize and re-partition the grids after importing. In this paper, we utilize
the said capability to an extent, importing a mesh and optimizing it using the
external MMG3D library (Dapogny et al., 2014). Considering the outputs of
YALES2 simulations, HDF5 format is used along with in-solver optimizations
which redistribute the input/output between threads to minimize the MPI
communications.

###  2.2 Governing Equations and Numerical Methods

In this Section, we will present a brief description of the governing
equations relevant to the work presented in this paper, followed directly by
specifications of the numerical methods used to solve them. We solve the
isothermal, incompressible two-phase fluid flow described by the Navier-Stokes
equations in the so-called single-fluid formulation (see e.g. Delhaye, 1973) :

$\frac{\partial\mathbf{u}}{\partial
t}+\nabla\cdot\mathbf{u}\otimes\mathbf{u}=\frac{1}{\rho}\left(\nabla\cdot\left(-p+\mu\mathbf{D}\right)+\sigma\kappa\mathbf{n}\delta_{s}\right)+\mathbf{f_{g}}$
(1)

with the continuity equation

$\nabla\cdot\mathbf{u}=0.$ (2)

In the above, $\mathbf{u},$ $p$ $\mu$ and $\rho$ stand, respectively, for
liquid velocity, pressure, viscosity and density. The interface is denoted by
$S,$ its normal vector as $\mathbf{n}$ and its curvature as $\kappa;$ the
delta distribution $\delta_{S}$ is centered on $S.$ In case of a discontinuity
of any of the material properties between liquid phases, or nonzero surface
tension coefficient $\sigma$, the pressure may also be discontinuous on the
interface $S$. This assumption is characteristic for the single-fluid
formulation, although it is assumed that velocity is continuous on $S.$ The
gravity (body) force is represented by $\mathbf{f_{b}}.$ The term $\mathbf{D}$
is the rate of deformation tensor

$\mathbf{D}=\frac{1}{2}\left(\nabla\mathbf{u}+\nabla^{T}\mathbf{u}\right).$

Note: If needed, we denote material properties with ”l” and ”g” indices
standing for liquid and gas; due to the single-fluid formulation these do note
appear in (1) as will be explained below (see e.g. (10).

Restriction of the curvature to the interface $S$ is done using the distance
function $\phi$, defined for each domain point $\mathbf{x}$ as minimum of its
distance from the interface:

$\phi(\mathbf{x},t)\colon=\min(d(\mathbf{x},S)),$ (3)

where $d$ is a chosen metric. We thus note
$\delta_{S}\colon=\delta_{S}(\phi)$. To denote restrictions to the interface
we may also use the subscript $S,$ such as
$\kappa_{S}\colon=\kappa\delta_{S}.$ When formulating the two-phase Navier-
Stokes equations (1), the singularity of the interface $S$ is solved thanks to
the the pressure jump, denoted using the jump notation
$[\cdot]_{S}=|\cdot_{l}-\cdot_{g}|.$ We thus define the pressure jump as

$[p]_{S}=(p_{l})_{S}-(p_{g})_{S}=\sigma\kappa_{S}+2[\mu]_{S}\mathbf{n}^{T}\cdot\nabla\mathbf{u}\cdot\mathbf{n},$
(4)

where $[\mu]$ is the viscosity jump if applicable. In all the above, the
normal vector $\mathbf{n}$ can also be defined using the distance function as

$\mathbf{n}=\frac{\nabla\phi}{|\nabla\phi|}.$ (5)

In the traditional level-set method formulation, the advection of the distance
function $\phi$ defined in (3) can be realized by a direct solution of its
advection equation

$\frac{D\phi}{Dt}=0$

, that is zeroing of the material derivative. This is made possible by the
fact that $\phi$ is continuous and smooth at least in the interface vicinity.
However, numerous works have indicated that the method is lacking in terms of
mass conservation of the traced liquid, especially in comparison to the
Volume-of-Fluid (VOF) type schemes (such comparisons can be found e.g. in
Aniszewski et al., 2014). Therefore, this work employs the Accurate
Conservative Level-Set (ACLS) method implemented in the YALES2 framework as
described in Janodet et al. (2022). In this method, a new scalar field $\psi$
is defined, using the hyperbolic tangent of the distance function:

$\psi(\mathbf{x},t)=\frac{1}{2}\left(1+\tanh\left(\frac{\phi(\mathbf{x},t)}{2\epsilon}\right)\right),$
(6)

where $\epsilon$ sets the profile thickness by modifying the steepness of
$\psi$ near the interface. By convention, $S$ is now defined by

$S(t)=\left\\{\mathbf{x}\in\mathbf{R}^{3}\colon\psi(\mathbf{x},t)=\frac{1}{2}\right\\},$
(7)

which is equivalent to setting $\phi=0$ in most cases333Equation (7) allows us
to count ACLS among the sharp interface methods, although it is interesting to
note that the definition of $S$ can be broadened by reusing the VOF convention
of considering as interfacial all grid cells for which
$0<\psi(\mathbf{x},t)<1$; this results in a diffuse interface method.. Due to
similarity to the $C$ function used in the VOF methods, we may occasionally
refer to $\psi$ as the volume fraction function although the definition (7) is
not identical to VOF volume fraction functions. However, $\psi$ assumes the
value $0$ in one phase and $1$ in the other, as $C$ does, which justifies this
naming.

Having defined $\psi,$ we can advect $S$ by a solving the equation:

$\frac{\partial\psi}{\partial t}+\nabla\cdot\left(\psi\mathbf{u}\right)=0.$
(8)

In YALES2, the solution of (8) is achieved via 4-th order spatial
discretization and the TFV4A time advancement scheme (Kraushaar, 2011). As is
the case with the LS method, ACLS requires re-initialization of the hyperbolic
tangent function, which in our work uses the approach of Chiodi and Desjardins
(2017). Like the normal vector $\mathbf{n}$, the curvature $\kappa$ is
calculated directly from $\phi$ (see Janodet et al., 2022 for the
implementation details).

Using $\psi,$ we are able to define the material properties for use in (1). We
utilize the definitions from the Volume-of-Fluid type methods (e.g. Xie et
al., 2017):

$\mu(\mathbf{x},t)=\mu_{g}+\left(\mu_{l}-\mu_{g}\right)\psi(\mathbf{x},t),$
(9)

where $\mu_{l}$ and $\mu_{g}$ are constant viscosities of liquid and gas,
respectively. With this smooth $\mu$ definition, intermediary velocity values
will be encountered in the near-interface grid nodes. Unlike (9), the
definition of $\rho$ is a sharp one and utilizes the Heaviside function
$H(x)$:

$\rho(\mathbf{x},t)=\rho_{g}+(\rho_{l}-\rho_{g})H\left(\psi(\mathbf{x},t)-\frac{1}{2}\right).$
(10)

This ensures that node-to-node density jumps will be coupled with the phase
definitions provided by (6), thus limiting the risk of non-physical momentum
fluxes (Vaudor et al., 2017).

$\star$

Considering the solutions of equations (1)-(10), we first note that YALES2 is
a projection method (Tryggvason et al., 2011) based solver. This means that
the Poisson equation for pressure will be solved; in this work this happens
via the Deflated Preconditioned Conjugate Gradient (DPCG) method (Malandain et
al., 2013) used to solve the resulting linear system. Inclusion of the
pressure jump (4) in the above solution is done via the Ghost Fluid Method
(GFM) which in our case is an unstructured grid formulation of Fedkiw et al.
(1999).

Similar to the coupled Level-set Volume-of-Fluid (CLSVOF) methods (Aniszewski
et al., 2014) the coexistence of the hyperbolic tangent function $\psi$ and
the distance function $\phi$ in ACLS leaves room for some choices in the
implementation. The functions can be advected independently and coupled (the
CLSVOF analog is presented in Aniszewski et al., 2012) or they can be re-
initialized periodically (or recreated) one from another (see Ménard et al.,
2007). In the YALES2 implementation of ACLS, the $\psi$ function is advected
by solving (8) and $\phi$ is actually recreated from $\psi$ at each solver
time-step. Then, the distance function $\phi$ can be used to calculate
$\mathbf{n}$ and $\kappa$. While the details of the $\phi$ calculation
algorithm can be found in Janodet et al. (2022), we summarize it briefly here.
First, YALES2 uses a narrow band of grid elements on both size of $S$ to
define $\phi$ only there, which reduces the CPU cost444The thickness of this
narrow band is a user-controlled parameter. Second, the reconstruction of
$\phi$ from $\psi$ is done using a geometrical Fast Marching Method in the
framework of the GPMM (Geometric Projection Marker Method) of Janodet et al.
(2019). It involves the full 3D triangulation of the interface on tetrahedral
grids, which is followed by projecting interface markers onto the triangulated
surface. The line along which each marker was projected is the first
approximation of the distance function gradient. Hence $\phi$ can be
reconstructed from this information by propagating the triangulated surface
via the Fast Marching Method inside the narrow band. This algorithm is
optimized for massively parallel calculations as described in Janodet et al.
(2022).

In closing of this Section, we denote that the distance function
$\phi(\mathbf{x},t)$ is essential for this work as it can be used to directly
obtain the distance between $S$ and selected grid nodes. For a point
$P(x_{p},y_{p},0)$ where $x_{p}$ and $y_{p}$ are coordinates on the initially
flat liquid interface, the elevation $\eta(t)|_{P}$ is approximately equal to
$\phi(P,t)$ which we use to generate the signal of the point probes (see Sect.
4.2). In general $\eta$ and $\phi$ can only be exactly equal if the local
curvature radius $R_{p}$ at $P$ will be greater than the distance between $P$
and $S$:

$\phi(P,t)=\eta(t)|_{P}\Leftrightarrow d(S,P)\ll R_{p}.$

While this condition is fulfilled for the probes used in this work, it has to
be noted that even the failure to observe it results in errors that are of the
order of the grid size near $S$, which is one to two orders of magnitude lower
than the $\eta$ signals measured.

###  2.3 Wave Generation

The wave generation technique used in this paper is based on the work (Rapp
and Melville, 1990) which introduced the multi-modal wave packet description
using the velocity potential $\Phi$. Unlike the solutions based on the piston
model (Chen et al., 2017) or momentum sources (Ha et al., 2013), this approach
is well suited for solitary or uni-modal waves. First, the explicit expression
for the liquid elevation writes as

$\eta(x=0,t)=\sum_{i-1}^{N}a_{i}\cos[\Omega_{i}(t-t_{f})-k_{i}(x-x_{f})],$
(11)

where $N$ is the number of modes. Formula (11) is not potential-dependent, and
can easily be implemented as a Dirichlet-type inlet condition for the
hyperbolic tangent function (6)555As a matter of fact, it is possible to
define naive wave generators using only (11) and a flat velocity profile,
however these simple techniques do not perform well in realistic loading
studies.. This is followed by velocity components defined as derivatives of
$\Phi$:

$u(x=0,t)=u_{0}+\frac{\partial\Phi}{\partial x}\quad
w(x=0,t)=\frac{\partial\Phi}{\partial z},$ (12)

while $v(x=0,t)=0.$ The system is closed by defining the potential

$\Phi(x,z)=\sum_{i=1}^{N}\frac{a_{i}r(t)\Omega_{i}}{k_{i}}\frac{\cosh(k_{n}(z+d))}{\cosh
k_{n}d}\sin[\Omega_{i}(t-t_{f})-k_{i}(x-x_{f})].$ (13)

In the above equations, $a_{i}$ $k_{i}$ and $\Omega_{i}$ are respectively the
amplitude, wavenumber and the angular frequency of each mode. Water depth is
represented by $d.$ The $t_{f}$ and $x_{f}$ characterize the focusing location
of the wave packet. The above equations are imposed as the time-dependent
Dirichlet boundary condition at $x=0.$ Additionally, the numerical ramping
function $r(t)$ is used in (13), defined as

$r(t)=2\arctan(8t)/\pi$

which ensures stability of the generator in the initial stage of the flow
(note $r(t)>0.96$ for $t>2$).

The system (11)-(13) has been proposed to produce multi-modal packets with a
set focusing location e.g. to generate breaking waves. However, as described
in Sect. 3, the simulated configuration consists of a series of identical
waves which interact with the cylinder cart. Therefore, we have decided on
setting $N=1$ in (12), and setting $a_{i}$ and $k_{i}$ to the values measured
experimentally. Additionally, we include $u_{0}$ in (12) corresponding to the
cart velocity (for the h5c3 configuration, $u=0.837$m/s). For $u_{0}\neq 0,$
using the ramping function $r(t)$ ensures the relaxation effect similar to
that of a fixed, near-inlet relaxation zone used e.g. by Qu et al. (2021),
except in the reference frame of the moving cylinder. It is easy to check that
assuming $N=1,$ equation (12) reduces to

$\displaystyle u(x,t)$ $\displaystyle=$ $\displaystyle
u_{0}-r(t)a_{1}\Omega_{1}\cos(kx-t\Omega_{1}-p_{s})\cdot$
$\displaystyle\cdot\frac{\cosh(k_{1}z+dk_{1})}{\cosh(dk_{1})}\mbox{ }$

and

$\displaystyle w(x,t)$ $\displaystyle=$
$\displaystyle-r(t)a_{1}\Omega_{1}\sin(k_{1}x-t\Omega_{1}-p_{s})\cdot$
$\displaystyle\cdot\frac{\sinh(k_{1}z+dk_{1})}{\cosh(dk_{1})},\mbox{ }$

where $p_{s}$ is the phase shift term equal to $\Omega_{1}t_{f}+k_{1}x_{f}.$
Since the wave generator matches the experimental parameters as explained in
Sect. 3, the wavenumbers and frequencies visible in (2.3) and (2.3) must obey
the dispersion relation (17) for the waves observed in the flume. With (17),
the above system is closed and can be implemented directly for any value of
$N.$

###  2.4 Post-processing

The YALES2 solver allows for the simulation data acquisition in multiple
formats (ASCII and binary), thus we will restrict this description solely to
the means used in this work.

The acquisition of the simulation scalar and vector fields is realized in
parallel using the HDF format. The solver uses HDF5 both for structured and
unstructured data; automatic balancing is used to distribute the I/O
operations between the used threads. However, the data chunks saved onto the
disks need not reflect the structure of the MPI threads. So for example, it is
perfectly possible to use 1024 processors which dump the data to 28 HDF5
chunks. Further processing of the data is possible using widely known
visualization suites such as Paraview (Ahrens et al., 2005). A priori all the
scalar and vector fields are periodically dumped during the simulation.

To offer a more fine-grained insight into the flow physics, we employ probes.
These are either of the predefined sampling points, lines or planes, over
which the acquisition is taking place at the specified time intervals. This
allows extraction of characteristics such as profiles or temporal evolution of
given scalar/vector fields (see e.g. Figure 6). In the simulations presented
here, we have opted for the point probes P1-P9 which we use to extract the
interface elevation $\eta$, local velocity magnitudes or pressure; details of
the placement of the probes are given in Table 2 (see also Section 4.2). Where
possible, the YALES2 probes are meant to correspond exactly or closely to the
placement of the resistive probes in the experimental setup, as explained
below. (When mentioning the experiment probes, we will denote them with the
”exp” subscript, such as P$3_{exp}$. )

Other information means are offered by several kinds of statistical
quantities. Starting from the lowest order statistics: the minimum, average
and maximum values of chosen scalar fields are stored as they evolve in time.
Typical uses here involve quantities such as $\max(\mathbf{u}(t))$ or
$(\int_{V}E_{K}dV)/V$. Moreover, YALES2 computes the mean (expectation) of
$\mathbf{u}$ understood as a random variable. Included are also higher moments
such as standard deviation (the root mean square, or the r.m.s.). This can be
performed also for other scalar fields such as the interface tracking
functions, and will be discussed in Section 4.7. Finally, the solver computes
statistics using the phase averaging procedure, as will be specified in
Section 4.7. Some statistical values, as well as quantities describing the
code performance, can be acquired from the solver logs.

For the majority of the data presented in this work, BASH and Gnuplot scripts
are used to automatize the presentation of the results, especially concerning
the ASCII data (probe signals, time traces). This post-processing involves
also the Fourier transforms (see Section 4.5) employing the FFTW3 (Frigo and
Johnson, 2005).

###  2.5 Mesh Generation

In this subsection, we will briefly describe the procedure of generating the
computational grid mesh(es) for the simulations presented in this work.
Additionally, a description of the actual used meshes will be given. The grid
generation procedure involves a specially constructed YALES2 test case whose
purpose is the import of a ASCII STL file describing the solid immersed in the
computational domain with arbitrary physical dimensions. In the procedure, an
ASCII file describing the solid is imported into YALES2 and placed in a domain
of arbitrary shape. Subsequently the STL geometry specified by the file is re-
meshed using tetrahedrons at a specified minimum spatial resolution. The user
controls parameters such as element skewness or the maximum mesh size gradient
$\max|\nabla(\Delta x)|$. In an intermediary stage we specify additional
refinement zones (.e.g. for the air-liquid interface), or coarse grid zones
(”sponges”) which may be used e.g. for numerical damping of kinetic energy and
testing purposes. Certain details of the used computational meshes are given
in Table 1.

Name | No. Elem. ($\cdot 10^{6}$) | $\min(\Delta x)$ |  $x$-span | Commentary
---|---|---|---|---
G0 | $21.8$ |  $5$mm |  $\pi$m | Tests, shorter cylinder, less refinement
G01 | $39$ |  $5$mm |  $\pi$m | Main production
G02 | $71.8$ |  $4$mm |  $\pi$m | Convergence study
G03 | $160$ |  $3$mm |  $\pi$m | Convergence study
G41 | $35.8$ |  $5$mm |  $\pi$m |  $C_{r}$ study, no cylinder
G42 | $39.7$ |  $5$mm |  $4.141$m |  $C_{r}$ study, longer, sponge layer
Table 1: Characterization of the computational grids (meshes) used in the
simulations. For all meshes, their physical dimensions are
$x\mbox{-span}\times 1\times 0.6$m. All grids are tetrahedral.  Figure 1: The
grid G01 as specified in Table 1. The grid is visible as a cut-volume along
the $xz$ plane. A - inlet, B - the fine mesh box containing the cylinder wake
area, C - the cylinder visible as a boundary fitted opening in the mesh. The
liquid/air interface location is presented with a dashed black line.

We denote that YALES2 employs the AMR (Adaptive Mesh Refinement) on run-time,
however in this work we use AMR only in the limited sense for grid generation
at the beginning of the simulations, the grid is then kept static.
Specifically, after the STL file describing the cylinder is imported and the
fine/coarse mesh zones prescribed, a preliminary simulation is ran which uses
AMR over several dozens of solver iterations until the mesh converges to
specifications. In this preliminary simulation, the MMG3D library (Dapogny et
al., 2014) is used to perform the grid adaptation (see also Benard et al.,
2016). Subsequently, thus obtained mesh is moved into a separate, production
simulation. The rationale for this procedure is that in our particular flow
case it is possible to predict e.g. the interface location in the quasi-
stationary flow stage. Even once the flow is fully developed and quasi-
stationary, it involves waves of known amplitude (approximately $2a$) thus the
interface resides in a zone of given thickness ($z$-span). We additionally
refine the vicinity of the solid cylinder and a box-shaped volume in the
cylinder wake, as shown in the Figure 1 which showcases the mesh G0 as
specified in Table 1. One notes by inspecting Figure 1 that any possible air
entrapment in the cylinder wake will take place in the refined zone as well.
This static-AMR approach circumvents the need of a fully fledged AMR on an
unstructured mesh thanks to the predictability of the interface residence
areas – this simplifies the simulation procedure and improves the
repeatability of results. Note that simulations of several other two-phase
flow cases have been published using YALES2 (Janodet et al., 2022) which
operate the AMR in a fully dynamic manner.

Of the grids described in Table 1, the grid G01 has been used as a production
grid providing majority of the results presented in this paper, except when
indicated. The G0 grid was used only in the testing stage, and was, like G01,
characterized by the $3.14$m domain length, however a different cylinder shape
was used than in G01, also the refinement areas in the cylinder wake were
significantly smaller, limiting the G0 grid to $21$ million elements. The
grids G02 and G03, carrying respectively $71.8$ and $160$ million elements,
are designed for the convergence study. Their respective minimal spatial
resolutions are $4$mm and $3$mm, otherwise they have identical specifications
to the G01 grid. Finally, two grids have been prepared as a part of the
investigations of the wave reflection phenomena (see Sect. 4.6) and
calculations of the reflection coefficient $C_{r}.$ The first of the grids,
G41 contains the channel geometry without the immersed solid cylinder, which
is a rudimentary method to ascertain the degree to which the latter
contributes to the (measured) wave reflections in the simulations. The
resolution is kept equal to the main production grid G01 and the number of
elements is $35.8\cdot 10^{6}.$ The other conceivable source of the
numerically-induced reflections is the outflow boundary condition, a
possibility we investigate using the last of the grids presented in Table 1:
G42. Compared to the main grid G01, the G42 has been elongated by an
additional $1$m downstream. Moreover, the additional domain volume thus
obtained acts as a numerical sponge layer, since it is left unrefined – which
explains a relatively low element count of G42 ($39.7$ million) compared to
G01, with only $7\cdot 10^{5}$ elements added. Adding sponge layers is a
rudimentary method of (numerically) absorbing the kinetic energy, and has been
well tested in wave hydrodynamics simulations (Ha et al., 2013; Israeli and
Orszag, 1981) including the wave-cylinder interactions (Ai et al., 2019). In
the case of the G42 grid, the unrefined elements in the sponge area have sizes
close to $\max(\Delta x)=0.05$m; no additional modifications of the velocity
field or the interface are used.

##  3 Experiment Specification

In this Section, we present a brief description of the experimental setup. The
experiments have been performed using a cylinder of $D=5$ cm made of
polymethyl methacrylate (PMMA) and clogged at its free-end. The cylinder is
mounted on a traction cart moving through a wave flume of 34 m in length,
$W=90$ cm in width and 120 cm in height filled with $d=0.4$ m of tap water.
The cylinder is partially immersed with a depth $h_{d}=0.23$ m in tap water,
and moves along the flume with a velocity, $u_{0}=0.8375$ m/s.

Figure 2: Configuration of the experiment with the cylinder of diameter $5cm$
and the resistive probe allowing the use of the point method of determination
of the adimensional coefficients.

Previously, a set of experiments, where cylinders have been translated at
different constant velocities through water initially at rest has been
achieved with the similar setup. For those experiments, the focus had been
brought on the drag and air-entrapment phenomena for different cylinder
diameters (Ageorges et al., 2019). In the current study, the waves are
generated using a piston type paddle. The incident/reflected wave amplitudes
and associated reflection coefficient $C_{r}$ are estimated using the method
of Mansard and Funke (Mansard and Funke, 1980) with two resistive probes
placed upstream the cylinder (see also Figure 2).

In the frame of reference of the laboratory, the incoming waves have the
period $T=2$ s, the amplitude (crest to through) $2a=0.048$ m and a maximum
celerity at free-surface $u_{FS}\simeq 0.127$ m/s (calculated from Stokes
theory at first order). The characteristics of the incoming waves are
determined in the frame of reference of the laboratory with probes placed at
P$3_{exp}(-1.5,-0.3,0)$ and P$2_{exp}(-0.7,-0.3,0)$ of the fixed cylinder and
a sampling frequency of $100$Hz. Note: these probe locations correspond to the
probe points specified in the Table 2 for the simulations.

In Figure 2, we present the configuration for the experimental study of the
problem. For the reminder of this description, we will, when specified, use
the frame of reference attached to the cylinder which is able to move with the
velocity of current $u_{0}$. A resistive probe is visible in both pictures
next to the cylinder at P1 position (see also Table 2). At the upper part of
the Figure 2, a piezoelectric sensor (Kistler 9327C), with a sensitivity of
-7.8 pC/N allows the measurements of the drag force $F_{x}$, acting on the
whole cylinder. The sensor is connected to a charge amplifier (Kistler 5015A)
equipped with a low pass filter of 30 Hz, itself linked to a LabView
acquisition card. Signals are filtered using a moving average on 20 points.
The sensor and the resistive probe are synchronised and have a frequency
acquisition of 100 Hz.

The relation between frequency and wave number is given by the dispersion
relation for capillary-gravity waves:

$\Omega^{2}=gk\tanh(kh),$ (16)

for $k<<k_{c}$ with $k=\frac{2\pi}{\lambda}$ the wave number, $k_{c}\simeq$
367 m-1, the capillary wave number, $\Omega=\frac{2\pi}{T}$ the pulsation, $g$
the gravity. In the frame of reference of the laboratory, the wave number is
equal to $k=1.7005$m-1.

In the reference frame attached to the cylinder adding a current by moving the
cylinder at constant velocity $u_{0}$ leads to the modification of the
dispersion relation (16). This Doppler effect modifies it via a frequency
shift of $ku_{0},$ namely

$\left(\Omega+ku_{0}\right)^{2}=gk\tanh(kh)$ (17)

leading to a period of 1.43 s in the cylinder frame of reference. The same
dispersion relation is implemented in the numerical wave generator, so that
the potential equation (12) can be closed.

##  4 Simulation Results

###  4.1 Specification of the Test Cases

The simulation utilizes a domain whose physical dimensions are meant to
closely resemble the actual experimental structure found at the LOMC. YALES2
makes it conceivable even to include the entire experimental geometry: the
$34$m-long tunnel, a moving measurement cart and a wave-generating beater –
since the solver includes provisions for moving solids with attached grids,
adaptive mesh refinement and is generally formulated along the LES methodology
(Pope, 1996) so that arbitrary precision levels can be chosen throughout the
domain. However, the associated CPU and GPU cost of such a holistic approach
would be prohibitive and unjustified, as there is hardly any need to e.g.
study the far-wake phenomena in the presented case. Thus, these authors have
opted for a more restricted study, which on the other hand allows for higher
spatial resolutions, the DNS-type methodology and inclusion of all transient
effects.

More specifically, the simulation domain is built using the frame of reference
of the cylinder: it remains at the point $(0,0,0)$ while the liquid flows
around it, with an inlet and outlet at both domain extremities. The cylinder
is placed along the $z$ axis, perpendicular to liquid interface at rest which
stretches in the $xy$ plane. The flow is from $x-$ to $x+$ while the $y-$ and
$y+$ domain walls are referred to as ”front” and ”back”. The sizes of the
domain are chosen as $3.14\times 1\times 0.6$m which corresponds to a
$3.14$m-long section of the experimental flume. The inlet is placed at
$x=-1$m, thus situating the outlet at $x=2.14$m. Water depth (zero-level, i.e.
disregarding the wave crest/through) is $d=0.4$m. Unless noted otherwise, the
results presented here correspond to the h5c3 experiment configuration, which
is a naming convention signifying the cylinder diameter $D=0.05$m, presence of
water current with the mean inlet velocity $u_{0}=0.8375$m/s, and the wave
amplitude (crest to through) $2a=0.048$m. For the h5c3 configuration, we
denote the values of the Keulegan-Carpenter number of

$K_{c}\colon=\frac{T\left(\mathbf{u}_{FS}\right)}{D}\approx 3.5,$

while the Reynolds number ${\rm Re}\,=41875,$ resulting in the ratio ${\rm
Re}\,/K_{c}\approx 11980.$

It is essential to note that the solid cylinder is, as in the experimental
setup, immersed only partially to the depth $h_{d}=0.23$m, thus leaving $17$cm
of spacing between itself and the channel (flume) bottom. Since our goal is to
simulate the moving measurement cart in its reference plane the boundary
conditions for the front and back walls are translating wall (translation
velocity set to $u_{0}$). No-slip walls are chosen for the bottom and top
walls as well as the cylinder surface. A free outflow is configured for the
$x+$ domain outlet, while the inlet condition is coupled with the wave
generator, as explained in the Section 2.3. A good overall look at the domain
shape is provided by Figure 5, in which the flow direction is from top to
bottom, and the cylinder (not rendered) is marked by an opening in the liquid
interface, the later coloured by $|u|.$

Considering the associated CPU cost of the presented simulations, it has been
controlled mostly by the CPU architectures/generations of choice. Using the
G01 grid ($39$ million elements) as reference , the cost ranged from $50\cdot
10^{3}$ CPUh when using Intel CPUs (Broadwell x86) down to $25\cdot 10^{3}$
CPUh using AMD (EPYC2 x86 Rome) cores – for a simulation spanning approx.
$20$s of physical time. Breaking down the algorithm cost per iteration results
in a even split between the LS advancement and redistancing ($\approx 45\%$
per iteration) and the Poisson equation solver ($\approx 43\%$). Looking in
turn at the cumulative (total) cost of the presented YALES2 simulations, the
i/o operations have a significant share in the CPU cost at approximately
$30\%$ (this however is strongly hardware-dependent).

###  4.2 Macroscopic Flow Evolution

Figure 3: Macroscopic flow geometry for the YALES2 simulation of cylinder-wave
interaction.

The macroscopic evolution of the interface for an example simulation is shown
in Figure 3. The surface is coloured according to the magnitude of the
velocity. The vertical cylinder is not rendered visible in Fig. 3, except as
an opening in the interface. The wake is however clearly visible behind it.
The cylinder position of the structure is $1$m from the inlet. In Figure 3, we
are observing an incoming wave upstream, the ”v-shape” liquid folds at the
borders of the cavity, as well as the traces of the vortical structures on the
liquid surfaces downstream. More accurate readings of the interface elevation
$\eta(\mathbf{x},t)$ are given below using the probes. At this stage, we only
denote that the result presented in Fig. 3, obtained using the $5$mm grid G01
closely resembles the macroscopic flow images obtained experimentally at this
Reynolds and Froude numbers (Ageorges, 2019). The level of detail e.g.
considering the interface geometry is far superior to certain previous works –
for example (Han et al., 2015) where the authors have used even finer meshes
near the immersed structure666The cited work concerned lower Reynolds numbers
in somewhat different context, but the obtained geometrical information
concerning the interface, as well as the details of the $\mathbf{u}$ field,
are much less detailed. and compares well to works such as Koo et al. (2014)
whose authors have investigated individual vortical structures around the
immersed structure at high Reynolds numbers. The same grid had been used for
the result present in Figure 4 below.

Figure 4:  The geometry of the near-cylinder interface for the same time
($t\approx 24.8$s) as in Fig. 3.

Figure 4 shows the same moment in the example simulation as in Figure 3,
except that this time we focus on the vicinity of the cylinder. The cylinder
is visible in blue, the interface is made semi-transparent, showcasing the
liquid elevation at the cylinder’s upstream side and a slight elevation drop
downstream. The plane visible in the background is coloured using the
vorticity norm or $|\omega|$, showcasing rich vortical structures in the wake.
The authors of Gonclaves et al. (2015) whose work concerns an experimental PIV
measurements for a case without waves777 However, it features a flow
characterized by ${\rm Re}\,\approx 43000$ which is close to our
configuration, even if the $h_{d}/D$ ratio is lower ($h_{D}/D\approx 2$
compared to $h_{D}/D\approx 4.6$ here. have established that the vorticity is
produced principally in two areas: one below the cylinder’s bottom, and
another present as a circulation zone directly adhering the cylinder
downstream. The action of the latter circulation zone is well visible in Fig.
4 due to high levels of vorticity.

Concerning $\omega,$ the authors of Ageorges et al. (2019) have established
that the air entrapment in the vicinity of the partially submerged structure
acts as to mitigate the vorticity generation in its wake; however their
findings were for ${\rm Re}\,>100000,$ while the case discussed here is barely
in the intermittent entrapment regime meaning that the bubbles are either not
at all, or rarely visible. Nonetheless, unlike their work, we have focused on
the wavy flow, which itself is a mechanism modifying $\omega$ as will be
discussed below in context of Figure 9. As said the macroscopic evolution of
the flow resembles closely the experimental results. The wake of the cylinder
evolves during the passage of the waves, we also observe the intermittent
entrapment of the air reappearing with the frequency $1/T$; in Figure 4 this
is visible in the form of small, elongated bubbles colored blue near the
cylinder.

Figure 5: Macroscopic flow geometry for the YALES2 simulation of cylinder-wave
interaction: top view, coloured by $|u|$. The flow direction is from top to
bottom.

Figure 5 presents the top-camera position view of the entire interface
geometry in the $3.14$m - long flume at the instant of $t\approx 3$s.
Colouring is done using the velocity magnitude. An incoming, generated wave is
clearly visible close to the inlet, while the cylinder wake is well visible
below it. The liquid folds forms a neck behind the structure (the ”v-shape”)
whose angle can be measured both for the experimental as numerical data (see
Figure 14). This instantaneous result presents also slight interface
deformations (visible at the outflow edge). This numerical phenomenon occurs
occasionally in YALES2 due to strong $\Delta x$ in the cells passed by the
interface (we noted it only in the G01 grid). This can be avoided using fully
dynamic AMR or finer meshes in general (such as G02 and G03), however even in
our meshing strategy these kind of errors appeared only penultimate cell layer
next to the outflow, and the error has no means to propagate upstream. Further
looks at the interface geometry will be provided below (e.g. Figures 12 and
15).

Figure 6: Positioning of point probes in the geometry of the YALES2
simulation. Each probe samples selected scalar fields every $n$-th time step
($\min(n)=1$). In this example, signals for the elevation of the $\eta$
interface are shown for four locations. In particular, the probe next to the
cylinder (number 1) is visible with its elevation signal showing $2a\approx
0.05$m, very close to the specification of the experimental configuration
h5c3..

Investigation of the water/air interface in more detail is made possible using
probes. For the current application, we have included a total of $9$ point
probes, all of which are placed at the $z=0$ plane. Probes measure any scalar
field chosen by the programmer; for this study we have settled upon $\eta$ and
$\mathbf{u}.$ Their positioning in the domain is portrayed in Figure 6, and
the full probe list with comments is provided by the Table 2.

Name | Position | Commentary | Experiment888Coordinates given in the cylinder frame of reference
---|---|---|---
P1 | $(-0.025,0.27,0)$ | cylinder proximity | P1${}_{exp}(-0.025,0.27,0)$
P2 | $(-0.7,-0.3,0)$ | near-inlet | P2${}_{exp}(-0.7,-0.3,0)$
P3 | $(-0.99,-0.3,0)$ | inlet | P3${}_{exp}(-1.5,-0.3.0)$
P4 | $(1,-0.3,0)$ | wake waves |
P5 | $(2.1,-0.3,0)$ | outlet |
P6 | $(-0.65,-0.3,0)$ | pair with P2 |
P7 | $(-0.6,-0.3,0)$ | pair with P2 |
P8 | $(1.1,-0.3,0)$ | pair with P4 |
P9 | $(0.0,0.4,0)$ | ”near wall” |
Table 2: List of the point probes in the flow with their spatial coordinates
and commentaries. For three of the probes, references are given to their
experimental counterparts. Note that for the probe P3, the simulation probe is
placed actually slightly closer to the cylinder ($1$m upstream as opposed to
$1.5$m in the experiment) which is due to the simulation domain extending only
$1$m upstream.

In Table 2, apart from specifying each probe’s location, we comment briefly on
its intended role; e.g. measurements of the wave elevation in the direct
vicinity of the cylinder for P1. One observes that the P3 (inlet) probe is
positioned in fact $1$cm downstream, which helps it convey more information
about the flow physics then it would at $x=0$ which would mean a raw copy of
the inlet boundary condition defined in Sect. 2.3. The cylinder proximity
probes P1 and P9 are placed both at the $x=0$ line to coincide with the
location of the solid in the channel. Probe P1 is positioned exactly in the
same position as the resistive probe in the experiment, as shown in the both
images of Figure 2. The other probe, P9 is positioned in the near wall region,
where it is possible to measure the elevation while avoiding the direct
contact with the ripples generated at the cylinder. In other words, the probe
P9 allows us to control the characteristics of the waves as seen by the solid
since, before the impact, the waves are effectively two-dimensional. Probes
such as P2 and P4 enable, respectively, a degree of control of the wave
amplitudes upstream (but not directly at the inlet) and in the cylinder wake.
In the meantime, probes P6, P7 and P8 are paired with P2 and P4 which enables
etc. the calculation of the reflection coefficient.

In the reference frame of the laboratory, the elevation $\eta$ seen by the
cylinder probe P1 for configuration h5c3 will simply be a near-exactly
sinusoidal curve with $2a=0.048$m and $T=1.43$s. This corresponds to the
bottom-left curve in Figure 6 which by comparison is characterized by slightly
higher wave amplitudes with $a\approx 0.032$m. We denote that our simple wave
generation method induces a small drop in amplitudes downstream, compared with
that of the inlet with $a$ dropping by $10$ to $20$ percent between probes P3
and P1. This is caused by the fact that while the incoming waves are governed
by (2.3) and (2.3), there is a possibility of backflows just downstream the
inlet, which interacts with the incoming wave. We believe that the resulting
wave amplitude at the structure is diminished by the said interaction, thus,
we have purposefully rescaled the injected $a$ to match the experiment at the
structure. Continuing the discussion of the $\eta$ evolution in the domain;
one observes that signals from the probe P4 downstream resemble P1 very
closely, accounting for the phase shift which is a known consequence of the
structure interaction (Sutherland and O’Donoghue, 1998). Finally, we denote an
amplified elevation signal at P5, suggesting that an additional treatment to
damp the wave kinetic energy might be necessary should this numerical method
be applied e.g. to longer domains and/or multiple structures.

###  4.3 Drag Coefficient and the Pressure Force Analysis

The pressure force exerted on the cylinder by the waves and current is
calculated in the simulations by integration of the water pressure at the
cylinder surface $\Gamma.$ To match the experimental measurement, we limit
this computation to the $x$-component of the force, by projecting it only onto
the first component of the vector $\mathbf{n}_{\Gamma}$ normal to the
boundary:

$F_{x}=-\int\limits_{\Gamma}\left(\mathbf{n}_{\Gamma}\right)_{x}pd\Gamma.$
(18)

One notices also that in the examined test case, the span-wise and vertical
pressure force components $F_{y}$ and $F_{z}$, while nonzero due e.g. to the
separation of turbulent structures – will not only be orders of magnitude
inferior to $F_{x},$ but also occasionally change sign depending on the actual
wave phase. Concerning the implementation, since the unstructured YALES2 code
uses tetrahedral meshing and boundary fitting, $\Gamma$ is merely a subset of
the general boundary object and no approximations are needed to ascertain its
position, unlike in certain structural mesh-based solutions (see e.g. Chen et
al., 2017).

Figure 7: The pressure force $F_{x}$ integrated on the surface of the
partially immersed cylinder. Data from the experimental measurement (circles)
and numerical simulations: grids G01 (line), G02 (dashed) and G03 (dotted).
Note that the time values ($t>150$s) correspond to the physical time of the
carriage traction in the experiment. The insets present, from left to right,
two signal peaks around $156.3$s and $159$s, respectively, and a minimum
around $161.1$s.

A typical evolution of the pressure force during $11$s of the simulated time
has been compared with the experimental measurement in Figure 7. The
experimental values come from a stage at which the measurement cart has
completed its acceleration, and the steady traction with the cylinder
impacting subsequent waves took place. Three numerical signals are visible in
Fig. 7, obtained using grids G01, G02 and G03. A very good correspondence is
observed on all grids, as far as the signal peaks are concerned. We observe a
slight disparity concerning the signal minima, as the experimental signal has
its minima at $3$N, which is over-estimated by the simulations at around $2$N.
To explain this slight disparity, we note first that $F_{x}$ is controlled
mostly by the height $h_{i}$ of the cylinder immersed in the water. It could
therefore be inferred that the simulated waves have a minimally lower through
(by which we mean $\min(\eta)$) than those in the experiment.

Figure 7 contains results obtained using grids G01, G02 and G03 as specified
in Table 1. In the upper plot presenting the full signal span, we observe that
both peaks and the minima decrease monotonically in value from G01 to G03.
However, using the inset images (bottom image row in Fig. 7) we find that the
$F_{x}$ grid dependency is slightly more complicated when examined in detail.
Namely, around the peak at $156.25$s (image bottom left), we observe that the
maximum peak value is obtained with G02 ($\min(\Delta x)=4$mm) and not G01.
Further on, around the peak at $159$s, where G01 yields the highest $F_{x}$
values, the G03 grid ($\min(\Delta x)=3$mm) offers a better prediction than
G02. Finally, the signal minimum at $161.1$s is best approximated with G02,
G03 and G01 in that order.

Figure 8: The pressure force $F_{x}$ signal obtained in the G01 simulation:
raw signal (small dots) and the $10$-point walking average (line).

The $F_{x}$ signal has a rather noisy character, as evidenced in Figure 8. In
the Figure, we have presented, for the G01 grid, the ten-point walking average
of the force signal (continuous line) along with the raw signal (dots). It is
clear that the variation in the raw signal is significant from iteration to
iteration. This may suggest, naturally, that the pressure distribution on the
immersed solid varies in a similar fashion; indeed the variations of the order
of $0.1-0.5$N seem consistent with the turbulence levels in the boundary layer
at the cylinder.

To conclude the discussion, we remark that the experimental $F_{x}$ signal as
shown in Figure 7 – especially the bottom-right inset – clearly has an
additional oscillation mode of approximately $8$Hz not seen in the numerical
results. After careful consideration and having re-run the simulations to no
avail, we have focused our attention on the experimental measurement system to
trace the source of this signal. Having performed a series of dry runs (only
cylinder cart traction without the liquid) and ping-tests, we have indeed
traced a relatively weak cylinder resonance at a frequency of $7.8$Hz (see
also Ageorges, 2019). This signal is superimposed onto the $F_{x}$ measurement
visible in Figure 7 and cannot be filtered out in post-processing (weak
$7.5$Hz and $8$Hz peaks are also visible in the FFT of the experimental
$F_{x}$ presented below in Figure 11). Other characteristic of the
experimental $F_{x}$ measurement is the application of a low-pass filtering
algorithm (cutoff frequency of $30$Hz) which in turn removes part of the
signal noise. For the future research, we remark that the solver allows
inclusion of the cylinder deformation with adaptive meshing which could help
resolve the mode present in the experiment.

###  4.4 Vorticity evolution

Figure 9: Superimposed curves of the $F_{x}$ and the vorticity integral
obtained in the G03 simulation. The vertical axis unit is Newtons, however it
so happens that $\int\omega dV$ oscillates around $4\cdot 1/s$, meaning that
no scaling is used in the plot and the $y$ tics are accurate also for the
vorticity unit.

The temporal evolution of the cylinder wake vorticity is presented in Figure 9
(dashed line). It should be precised that in the numerical simulation results
we have found that the temporal evolution of global vorticity in the whole
domain is indistinguishable from vorticity integrals calculated over wake sub-
domain. In other words, the bulk vorticity is contained in the wake, thus we
are presenting the global $\omega$ integral in Figure 9. In order to be able
to trace the interaction between the major waves, the cylinder and the wake
vorticity, we have superimposed onto the plot the $F_{x}$ curve obtained from
the G03 simulation (dotted line in Fig. 7). As can be readily observed in Fig.
9, there is a strict synchronisation between the $\omega$ maxima and that of
the pressure force. Namely, two $\omega$ peaks are found for each wave period.
Thus, vorticity minima are observed whenever wave crests or throughs impact
the cylinder, with $\omega$ peaking whenever a wave slope is directly upstream
or downstream. This behaviour has previously been described e.g. by Arabi et
al. (Arabi et al., 2019) concerning an interaction of waves with a cubical
obstacle which they studied numerically. Upon the impact of the incident wave
at the windward cylinder side, vortical structures are formed there due to the
creation of the water elevation (damming). Also, at this stage the structures
near the cylinder bottom intensify. This leads to the creation of the first
$\omega$ peak, while the second peak forms after each incident wave’s crest
passes the cylinder, leading to the intensified creation of the detached
vortical structures on the leeward side. The vorticity evolution is affected
also by the reshaping of the cavity (change in the wake angle $\alpha$, see
Figure 14).

In context of the vortical structures detaching from the cylinder, we denote
that the velocity gradient in the boundary layers has a certain grid
dependency. This poses a challenge on homogeneous tetrahedral grids (see Mukha
et al., 2021). In the case of the G0-G03 grids used here, we have denoted slow
convergence of the wall-normal velocity gradient with the values of $Y+\leq
20$ near all boundaries (notably, $Y+<20$ next to the $z-$ wall (flume bottom)
as well as cylinder surface). However, both purely geometrical aspects (such
as the grid skewness) as well as flow physics (boundary layer separation,
cusps in geometry), may give rise to occasionally higher $Y+$ values. We have
denoted the $\max(Y+)\leq 90$ for all G01, G02 and G03 grids, with the maximum
values corresponding either to highly skewed tetrahedra, or the cusps e.g. at
the bottom of the immersed cylinder.

###  4.5 Spectral Analysis

In this section, we briefly present the results of the investigations using
spectral information obtained mainly by transforming the probe signal. In two-
phase flows, while the Fourier transforms (FFT, DFT) of the signals sampled
far from the interface are a well accepted tool (Pope, 1996), the number of
works employing it to actually study the interface itself is far inferior.
Ling et al. (Ling et al., 2015) have used it on interface height samples
obtained from a shearing flow. Our application is similar - for both cases the
key fact is that the interface is initially flat, so that the height (or
distance) - measuring probe can yield meaningful information on interface
wrinkles, oscillations or deformations as it always oscillates around the
initial zero-height position. In applications such as atomization (Ling et
al., 2015) or even Rayleigh-Taylor instability’s later stages where the bulk
of the liquid undergoes a complete displacement (Aniszewski et al., 2014) this
kind of application would not be possible.

Figure 10: Spectra of the $\eta$ signal for six of the probes using G03 (a-d)
and G41 (e-f) grids.

That said, Figure 10 presents the discrete FT on six elevation $\eta$ signals
obtained by the probes throughout the simulation. The spectra only supply
information for $f<10$Hz, as above that frequency, for this flow
configurations, no signals of amplitude stronger than noise are expected to
appear in the flow. In Figure 10, the spectra of the G03 simulation
(Fig.10a-d) and the G41 grid (no cylinder, Fig.10e-f) are visible. All spectra
are characterized by the dominant first harmonic ($T=1.43$s or $f=0.69$Hz, the
generator frequency). Unlike upstream (Fig.10b), all spectra downstream the
cylinder starting with the structure probe itself (Fig.10a,c,d) contain the
second harmonic near $2f=1.38$Hz. We could associate this frequency with the
peaks of $\omega$ as described in the context of Figure 9. The third harmonic
($3f=2.07$Hz) is well visible further downstream, suggesting it may originate
e.g. from the action of the evolving wake edges. Finally, the outflow probe
spectrum (Fig. 10d) has evidence of further harmonics. In contrast, the flow
for which there was no fluid-structure interaction (Fig. 10e-f) save for the
containment effect of the channel (Potgieter et al., 2022), is characterized
by lower spectral content limited to the third ($1$m downstream, Fig.10e)
and/or fourth (outlet) harmonics. Higher frequency content is largely absent
in this flow.

Figure 11: Top: Elevation $\eta$ spectrum of the P1exp-P3exp probes and drag
force $F_{x}$ (black line). Bottom: probe P1exp elevation (blue) and drag
force spectra (black) for the cylinder translating at the velocity
$u_{0}=0.8375$ m/s

The experimental spectra of the incoming waves (in the frame reference
attached to the laboratory), the drag force and the probe P1 attached to the
moving cylinder are depicted in Figure 11. In comparing the experimental
signal spectra presented in Fig.11 with Fig.10, we use the blue curve
presented in the upper picture of Fig.11, i.e. the $\eta$ signal associated
with the P1exp probe. This curve corresponds directly to Fig. 10a. Taking into
account the difference in the amplitude normalization, we notice a similarity
in the slope of both spectra. Also, the second harmonic visible in Fig.10a has
a far lower amplitude compared to Fig.11. Obviously, the experimental data
contains at least ten clearly visible harmonic peaks, and has excellent peak
to noise ratio, the lack of which presents an obstacle in interpreting
Fig.10a.

###  4.6 Analysis of Wave Reflection

Both numerical and experimental data in our study allow for estimations of the
characteristics of the reflection coefficient. In the simulation, probe
elevation signals are used, allowing for a free choice of the involved probe
pairs or triplets. Subsequently, we employ the methodology of Goda and Suzuki
(1976) to resolve the composite waves into those incident and reflected. The
calculation is performed in the spectral space, with the restriction that the
elevation records must originate from probes spaced closer than $\lambda/2.$
For the experimental results measured using the resistive probes (Mansard and
Funke, 1980), the exact same methodology to calculate $C_{r}$ is used.

Thus, for the experimental data in the h5c3 case, the reflection coefficient
of the incoming waves is estimated at $C_{r}\simeq 0.0451$ using probes P2exp
and P3${}_{exp}.$ In the simulations, the base grid G01 yields $C_{r}\simeq
0.35$ using probes P2 and P6. However, using the G03 with the same probes
involved, we have obtained $C_{r}\simeq 0.21,$ which is still roughly $5$
times the experimental value. It seems an expected result considering the
simplicity of our wave generation scheme and the lack of more rigorous outflow
damping. It has to denoted that the probe signal is – both for the
experimental and numerical data – truncated to remove first seconds of the
flow corresponding to the cart acceleration stage999Certain sensitivity to the
exact choice of the truncation time takes place, thus we report averaged
$C_{r}$ values here..

Additionally, we have performed studies using the grids G41 and G42. Results
using the G41 grid indicate that $C_{r}$ actually grew to $0.5$, confirming
that the reflection is produced at the outlet in all tested simulations (and
not only the default configurations with the cylinder G01-G03). Since the G01
simulation contains the cylinder and G41 does not, and both the measurement
points lie upstream the cylinder, we conclude that the reflected wave was
partly dispersed by the cylinder in G01-G03 before reaching the location of
the probes. This explains the growth of $C_{r}$ using the G41. Additional
simulations were performed using the G42 grid (see Table 1) which contains a
longer flume and applies a 1m long sponge layer, i.e. a coarsen grid with
larger tetrahedra. Unfortunately, the results for this numerical setup yield
$C_{r}>1$; this amplification can be explained twofold: first, the sponge
layer is insufficient to damp the outgoing waves using only the phenomenon of
numerical viscosity (i.e. numerically induced $E_{k}$ dissipation). Second,
the reflected wave moves back upstream in the G42, the action no more that
efficiently dispersed by the obstacle as in the case of G01-G03.

Thus, we conclude that obtaining an exact match of the $C_{r}$ with the
experimental data could be obtained either by additional damping techniques
near the outflow (e.g. porous ramps as in the laboratory setup), or the
further increase of the spatial resolution, in hope that the dependency of
$C_{r}$ on the mesh size – as observed progressing from G01 to G03 – continues
to the finer meshes.

###  4.7 Statistical Analysis

Instantaneous results limit our possibilities in discerning between transient
and periodic phenomena especially at the small scales. Thus, e.g. to further
analyze the possible presence of reflected components in the composite waves
present in the domain, we have employed statistical techniques such as phase
averaging. The technique is applicable to periodic processes (see Wernert and
Favier, 1999 or for a similar application Corvaro et al., 2019) and utilizes
temporally averaged data grouped in bins spaced by $T$ (in which case $T$ is
the inherent process period, which in our case is naturally the wave period).
More specifically for the time-varying signal such as the elevation $\eta(t)$,
the phase average can be estimated as

$\langle\eta(t,N)\rangle_{est}=\frac{1}{N}\sum\limits^{N}_{i=1}\eta(t+(i-1)T),$
(19)

where $N$ is the number of the averaged periods which should be as large as
possible. As explained by Wernert and Favier (1999), the desirable $N$ values
are of order of hundreds for the experimental signals. However, in this study
we are limited in terms of the CPU cost, limiting the temporal span of a
typical G01 simulation to $t_{end}\in[13T,30T]$ with $t_{end}$ being even
smaller for the more costly grids such as G03. Thus, results presented below
should be treated as rough approximations in terms of the convergence of (19)
to $\langle\eta(t,N)\rangle.$

Figure 12: The phase-averaged mean values of the volume fraction function
$\chi(\phi)$ for the simulation on the G03 grid. (The $\chi(\phi)=0.5$
isosurfaces are shown.) Each picture corresponds to a $T/10$-sized collection
bin, where $T=1.43$ (e.g. the top left image is obtained using data from the
bin $b_{1}$ where $t\in[0+nT,0.143+nT]$). The data has been collected over the
entire length of the G03 simulation (approx. $12$s of physical time, thus
$n=1\dots 8$).

An example of the phase-averaged mean estimation of the elevation is presented
visually in Figure 12, which showcases the isosurfaces of the
$\langle\phi\rangle_{est}$ for $N=8$ obtained in the G03 simulation. Since the
images in Fig. 12 correspond to all ten bins $b_{1}\dots b_{10}$ spanning the
entire period $T,$ we observe how the $\lambda\approx 2.34$m wave101010The
wavelength is given in the frame of reference associated with the cylinder. In
the laboratory frame of reference $\lambda\approx 3.69$. passes through the
cylinder proximity in the $b_{1}$ image, leaves the simulation domain in the
$b_{5}$ and $b_{6}$ bin images to reenter by the $b_{7}.$ Other notable
features include fluctuations in the cylinder trace (v-shape) angle, which we
will discuss below.

As made clear by Figure 12, the numerical setup results in a low-amplitude,
“w-shaped” reflection wave close to the outlet which remains nearly static
during the flow period $T.$ This is consistent with the relatively high values
of the $C_{r}$ as described in Sect. 4.6; we attribute this to the rather
short computational domain, and the lack of numerical wave suppression towards
the channel outlet. However, we were not able to confirm the possible
propagation of the reflected wave towards the immersed cylinder. Of similar
(numerical) origins are the two-dimensional, stripe-shaped disturbances
visible upstream of the cylinder in Figure 12. These disurbances result from
the small-scale disturbances at the domain inlet, and are more prominent when
$a$ (the injected wave amplitude) is higher. In most cases, however, as
visible in Figure 12 (see e.g. $b_{6}$) these formations dissipate before
impacting the structure, and their effect on data such as the $F_{x}$ spectra
presented in Sect. 4.5 is negligible.

Figure 13: The phase-averaged mean value of $\chi(\phi)$ for the G03
simulation. Image of the $\chi=0.5$ iso-surface presented previously in Figure
12, with the color scheme corresponding to the phase-averaged standard
deviation of $\phi$.

The isosurface corresponding to the third bin $b_{3}$ of the phase average
(Figure 12) is reprinted in Figure 13 with a different color scheme. This
time, the color corresponds to the root mean square
$\langle\phi^{2}\rangle^{1/2}$ of $\phi$ obtained with phase
averaging111111Referred to also as the standard deviation, i.e. the square
root of the variance, see Pope (1996), p.41.. Using such color scheme allows
us to focus our attention on the deviation of the interface position from its
phase-averaged mean position for the chosen bins (i.e. selected intervals of
the wave period). Figure 13 confirms, in this way, that the areas with lowest
period-to-period variance correlate with the crests and throughs of the
$\lambda=2.34$m waves. Interestingly, an area of the v-shape in the downstream
trace is also coloured blue in Fig. 13 suggesting that the v-shape evolution
is strictly periodic and not subject e.g. to turbulence in the cylinder wake
at least at the tested Reynolds and Froude numbers. It is also interesting to
note that the near-outflow corresponds to maxima of
$\langle\phi^{2}\rangle^{1/2}$, i.e. the variance of the interface position is
high near the outflow during the simulation.

Figure 14: The ”v-shape” half-angle study: experimental result (red squares)
and simulation (black triangles). The inset presents the angle-measurement
technique as explained in the text.

In order to further validate the numerical result using the phase averaging,
we finalize the section with a comparison between the numerical and
experimental results concerning the half-angle of the "v-shape", which is the
geometrical trace created by the liquid folds behind the cylinder. The
comparison is presented in Figure 14. The comparison is performed as follows.
For each of the bins, the $\langle\chi\rangle=0.5$ iso-surface, representing
the mean interface position, is examined. Using the top view (i.e. along the
$z$ axis), we first process the image with edge detection algorithms and
subsequently measure the angle between the domain $x$ axis and the liquid
folds (measurement is taken $2D$ downstream). The procedure is presented
symbolically in the inset of Fig. 14: we can see two diagonals marking the
liquid folds. (In this particular example (bin $10$), the folds are not
strictly symmetrical, in which case an average value is presented in Fig. 14.)
The error bars in the numerical case are obtained using the standard deviation
information for the $\chi(\phi)$ function.

Analyzing the Figure 14, we observe a relatively close fit between the half-
angle values for $7$ of $10$ bins. Notably, the wide ”v-shape” half-angle
$\alpha>35$ for the phase $4\pi/5$ is well captured by the simulation. This
moment in the phase-average corresponds to the onset of the impact of the
$\lambda=2.34$m wave on the cylinder, the actual wave crest passing through
the structure takes place for the next measurement point ($\pi$) which is also
well captured. The points $6\pi/5$ and $7\pi/5$, corresponding to the moments
in the flow where the large wave crest has moved downstream the cylinder,
exhibit discrepancies between simulation and experiment. As shown in the
vorticity analysis (see Figure 9) once the crest clears the cylinder a
vorticity peak appears, suggesting a maximum vorticity production in the wake.
Since this mechanism is severely under-resolved in our simulation, we suppose
that the vorticity production mechanism underwater may imprint on the surface
eddies in the wake influencing the half-angles. Thus, the discrepancies in
Fig.14 could be of numerical origin and due to lack of resolution. That said,
the general $\alpha$ evolution seems well predicted by the simulation within
these limitations, as evidenced by Fig. 14.

The work of Rabaud and Moisy (2013) proposes a method to analytically
approximate the (half) wake angles based on the Froude number. They show that
for ${\rm Fr}\,>{\rm Fr}\,_{c}\approx 0.49,$ the angle expresses as:

$\alpha=\tan^{-1}\frac{\sqrt{2\pi{\rm Fr}\,^{2}-1}}{4\pi{\rm Fr}\,^{2}-1},$
(20)

which is applicable to our configuration since ${\rm Fr}\,=1.2$, leading to
the ”average” wake angle of $\alpha=18\degree.$ Applying (20) to instantaneous
velocities probed in the interface vicinity with the $P_{1}$ probes, we
denote, that $u\in[0.4,1.6]$ there resulting in ”local” values of the Froude
number between $0.57$ and $2.28$, hence $\alpha_{min}\approx 14\degree$ and
$\alpha_{max}\approx 36\degree,$ which is in good comparison with Figure 14.

Figure 15: A visual comparison between instantaneous results of the
experimental (left) and numerical (ritght) study. The image corresponds to
$t\approx 152.6$s in the laboratory frame of reference.

Final comparison between the numerical and experimental results is found in
Figure 15, showcasing a rendered view of the simulation (right) with the top-
view photography obtained during the experimental acquisition (left). The
laboratory view has been obtained with a portable camera fitted directly to
the cylinder mount, with a similar (although not identical121212This is
evident from the fact that far more liquid is visible upstream the cylinder in
the rendered image.) position used for the rendering. The view corresponds to
an instantaneous interface geometry at $t\approx 152.6$s (in the laboratory
frame of reference), and a wider v-shape angle. In the photograph visible in
Fig. 15, the horizontal yellow-green lines correspond to the physical flume
walls, while in the rendering they are artificially superimposed (the actual
camera angle did not include domain walls at the $y$ extremities). Concerning
the interface characteristics, while the v-shape is comparable in angle, there
seems to be more of a structure visible in the wake $4-5$D downstream in the
experimental image. This may suggest that the sub-surface wake eddies are
similarly weaker in the numerical result, as noted in the context of Fig. 14.
The upstream liquid elevation appears in both images; also a similar interface
structure (small-length waves) is found in both images further downstream.

##  5 Conclusions

In this paper, we have presented a numerical tool for simulating the
interactions of the waves with a submerged structure. Simulations involve a
partially submerged cylinder towed in a flume containing wavy flow. Our study
involves an experiment wherein that exact setup is used in a full length
flume. The experimental campaign conducted at the LOMC laboratory includes
measures of multiple cylinder radii at varying ${\rm Re}\,$ and $Fr$ values;
in our paper the focus has been on a single configuration where $D=0.05$m and
${\rm Re}\,\approx 41875.$ The simulation domain contains only a section of
the flume, and the waves are generated at the inlet.

We have observed a very good simulation-to-experiment comparison concerning
wave geometry, macroscopic flow statistics and, most importantly, the
hydrodynamic load onto the cylinder expressed in terms of the pressure force
$F_{x}.$ From a numerical point of view, the faithfulness of the $F_{x}$
calculation seems to rely mostly on the correct immersion depth for the
cylinder, as well as the mesh quality. This yields the proper integration
surface on the solid boundary and improves the match with the experimental
data. This, subsequently, means that factors such as fitted boundary
smoothness (lack of any artificial roughness) and/or the quality of the
boundary tetrahedra may be equally influential as the minimum grid size
$min(\Delta x).$

We have found that the vorticity generation mechanism is closely linked with
the position of the incident waves, which is consistent e.g. with previous
studies (Ageorges et al., 2019, 2021). Each passing wave has two associated
$\omega$ peaks in the cylinder wake zone as described by Arabi et al. (2019) –
one at the moment of the impingement, and another one associated with the
vortical structure separation and the growth of the ”v-shape” angle $\alpha$ –
after it passes the solid structure. Our study can be considered highly
resolved in comparison to certain previous works (Corvaro et al., 2019).

Using the information from the probe points placed in the calculation domain,
we were able to establish a close correspondence with the experimental results
in the spectral space. A singular, not-matching peak was found in the
experimental data caused by the resonance of the cylinder traction cart. Data
from the flow probes has allowed to find a closely matched $\eta$ evolution,
as well as control phenomena such as wave reflections. Due to significant
values of the reflection coefficient $C_{r}$ measured by the probes close to
the inlet, we must conclude that using the $3.14$m-long domain results in
significant levels of reflection. A ’w-shaped’ reflected wave is visible in
the phase-averaged geometry plots, and its removal should be the natural task
for the continued research. In fact, during the work described in this paper,
longer computational domains (up to $11$m) have been tested, including
variants with outflow ramps. However, this solution was deemed too
computationally expensive for the current campaign. Separately, wave
absorption zones (Han et al., 2015) were tested, however the latter solution
did not mitigate the reflection phenomenon in preliminary tests. Thus, in
these authors’ opinion, a solution combining both a ramp – with a form of
simulated porosity – and absorption zones should be tested in a separate work.
The same goes for techniques such as fully dynamic AMR.

That said, the goal of the paper was to demonstrate that even with a
relatively short flume section simulated, it is possible to capture
fundamental wake physics, especially if the spatial resolution is significant
close to the embedded structure. Reusing the presented method can present an
interesting opportunity e.g. for testing higher Reynolds numbers, which in
turn would place the flow in the regime characterized e.g. by significant air
entrapment in the cylinder wake. In broader terms, it allows a more wide-
ranging parametric studies of the wave-structure interaction in a two-phase
flow.

## Acknowledgements

The authors acknowledge the funding in the frame of the LabEx EMC3-2019-PC
STRUCTIMM project at the CNRS. Calculations were performed using the resources
of the CRIANN (Centre Régional Informatique et d’Applications Numériques de
Normandie) notably the Myria supercomputer (2020-2021). More calculations were
performed using the TGCC (Très Grand Centre de calcul du CEA) supercomputer
Irene Joliot-Curie (2021).

## References

  * Ageorges (2019) Ageorges, V., 2019. Ecoulement et entrainement d’air autour d’un cylindre vertical partiellement immerge. Ph.D. thesis, Normandie Université, Le Havre, France.
  * Ageorges et al. (2019) Ageorges, V., Peixinho, J., Perret, G., 2019. Flow and air-entrainment around partially submerged vertical cylinders. Physical Review Fluids 4, 064801.
  * Ageorges et al. (2021) Ageorges, V., Peixinho, J., Perret, G., Lartigue, G., Moureau, V., 2021. Experiments and simulations of free-surface flow behind a finite height rigid vertical cylinder. Fluids 6 (10).
  * Ahrens et al. (2005) Ahrens, J., Geveci, B., Law, C., 2005. 36 - paraview: An end-user tool for large-data visualization. In: Hansen, C. D., Johnson, C. R. (Eds.), Visualization Handbook. Butterworth-Heinemann, Burlington, pp. 717–731.
URL https://www.sciencedirect.com/science/article/pii/B9780123875822500381

  * Ai et al. (2019) Ai, C., Ma, Y., Yuan, C., Dong, G., 2019. A 3d non-hydrostatic model for wave interactions with structures using immersed boundary method. Computers and Fluids 186, 24–37.
  * Aniszewski et al. (2012) Aniszewski, W., Boguslawski, A., Marek, M., Tyliszczak, A., 2012. A new approach to sub-grid surface tension for LES of two-phase flows. Journal of Computational Physics 231, 7368–7397.
  * Aniszewski et al. (2014) Aniszewski, W., Ménard, T., Marek, M., 2014. Volume of fluid (VOF) type advection methods in two-phase flow: A comparative study. Computers & Fluids 97 (0), 52 – 73.
  * Arabi et al. (2019) Arabi, M. G., Sogut, D. V., Khosronejad, A., Yalciner, A. C., Farhadzadeh, A., 2019\. A numerical and experimental study of local hydrodynamics due to interactions between a solitary wave and an impervious structure. Coastal Engineering 147, 43–62.
URL https://www.sciencedirect.com/science/article/pii/S0378383918302333

  * Benard et al. (2016) Benard, P., Balarac, G., Moureau, V., Dobrzynski, C., Lartigue, G., D’Angelo, Y., 2016. Mesh adaptation for large-eddy simulations in complex geometries. International Journal for Numerical Methods in Fluids 81 (12), 719–740.
URL https://onlinelibrary.wiley.com/doi/abs/10.1002/fld.4204

  * Benard et al. (2018) Benard, P., Viré, A., Moureau, V., Lartigue, G., Beaudet, L., Deglaire, P., Bricteux, L., 2018. Large-eddy simulation of wind turbines wakes including geometrical effects. Computers and Fluids 173, 133–139.
URL https://www.sciencedirect.com/science/article/pii/S0045793018301154

  * Bernard et al. (2020) Bernard, M., Lartigue, G., Balarac, G., Moureau, V., Puigt, G., 2020. A framework to perform high-order deconvolution for finite-volume method on simplicial meshes. International Journal for Numerical Methods in Fluids 92 (11), 1551–1583.
URL https://onlinelibrary.wiley.com/doi/abs/10.1002/fld.4839

  * Brocklehurst et al. (2011) Brocklehurst, P., Korobkin, A., Pǎrǎu, E. I., 2011. Hydroelastic wave diffraction by a vertical cylinder. Philosophical Transactions: Mathematical, Physical and Engineering Sciences 369 (1947), 2832–2851.
URL http://www.jstor.org/stable/23035834

  * Chen et al. (2017) Chen, Q., Zang, J., Kelly, D. M., Dimakopoulos, A. S., 2017. A 3-d numerical study of solitary wave interaction with vertical cylinders using a parallelised particle-in-cell solver. Journal of Hydrodynamics 29, 790–799.
  * Chiodi and Desjardins (2017) Chiodi, R., Desjardins, O., 2017. A reformulation of the conservative level set reinitialization equation for accurate and robust simulation of complex multiphase flows. Journal of Computational Physics 343, 186–200.
URL https://www.sciencedirect.com/science/article/pii/S0021999117303327

  * Cleary and Rudman (2009) Cleary, P. W., Rudman, M., 2009. Extreme wave interaction with a floating oil rig: prediction using sph. Progress in Computational Fluid Dynamics, an International Journal 9, 332–344.
  * Corvaro et al. (2019) Corvaro, S., Crivellini, A., Marini, F., Cimarelli, A., Capitanelli, L., Mancinelli, A., 2019. Experimental and numerical analysis of the hydrodynamics around a vertical cylinder in waves. Journal of Marine Science and Engineering 7, 1–18.
  * Dapogny et al. (2014) Dapogny, C., Dobrzynski, C., Frey, P., Apr. 2014. Three-dimensional adaptive domain remeshing, implicit domain meshing, and applications to free and moving boundary problems. Journal of Computational Physics.
URL https://hal.sorbonne-universite.fr/hal-00804636

  * Davis (1974) Davis, A., 1974. Short surface waves in the presence of a submerged circular cylinder. SIAM Journal of Applied Mathematics, 479–491.
  * Delhaye (1973) Delhaye, J., 1973. Jump conditions and etropy sources in two-phase systems. local instant formulation. International Journal of Multiphase Flow 1, 395–409.
  * Dubief et al. (2020) Dubief, Y., Balarac, G., Mendez, S., Mercier, R., Bégou, P., Benard, P., Lartigue, G., Moureau, V., 2020. High-fidelity simulation as a decision-making tool for the protection against aerosols carrying covid-19. In: Séminaire SMI de l’Académie des Science.
URL http://www.legi.grenoble-inp.fr/web/spip.php?article1599&lang=fr

  * Elhanafi et al. (2021) Elhanafi, A., Duffy, J., Macfarlane, G., Binns, J., Keough, S. J., 2021. Inline forces and bow wave height on a vertical cylinder moving in waves — experimental study and cfd validation. Journal of Fluids and Structures 107, 103387\.
URL https://www.sciencedirect.com/science/article/pii/S0889974621001705

  * Fedkiw et al. (1999) Fedkiw, R., Aslam, T., Merriman, B., Osher, S., 1999. A non-oscillatory eulerian approach to interfaces in multimaterial flows (the ghost fluid method). Journal of Computational Physics 152, 457–492.
  * Fondelli et al. (2015) Fondelli, T., Andreini, A., Facchini, B., 2015. Numerical simulation of dam-break problem using an adaptive meshing approach. Energy Procedia 82, 309–315.
  * Frigo and Johnson (2005) Frigo, M., Johnson, S. G., 2005. The design and implementation of FFTW3. Proceedings of the IEEE 93 (2), 216–231, special issue on “Program Generation, Optimization, and Platform Adaptation”.
  * Goda and Suzuki (1976) Goda, Y., Suzuki, Y., 1976. Estimation of Incident and Reflected Waves in Random Wave Experiments. pp. 828–845.
URL https://ascelibrary.org/doi/abs/10.1061/9780872620834.048

  * Gonclaves et al. (2015) Gonclaves, R., Franzini, G., Rosetti, G., Meneghini, J., Fujarra, A., 2015. Flow around circular cylinders with very low aspect ratio. J. Fluids Struct. (54), 122–141.
  * Ha et al. (2013) Ha, T., Lin, P., Cho, Y.-S., 2013. Generation of 3d regular and irregular waves using navier-stokes equations with an internal wave maker. Coastal Engineering 76, 55–67.
  * Han et al. (2015) Han, Y., Zhan, J.-M., Su, W., Li, Y., Zhou, Q., 2015. Numerical simulation of in-line response of a vertical cylinder in regular waves. Journal of Coastal Research 31 (4), 879–891.
  * Hu et al. (2016) Hu, Z. Z., Greaves, D., Raby, A., 2016. Numerical wave tank study of extreme waves and wave-structure interaction using openfoam®. Ocean Engineering 126, 329–342.
URL https://www.sciencedirect.com/science/article/pii/S0029801816303936

  * Israeli and Orszag (1981) Israeli, M., Orszag, S., 1981. Approximation of radiation boundary conditions. Journal of Computational Physics 41, 115–131.
  * Janodet et al. (2022) Janodet, R., Guillamón, C., Moureau, V., Mercier, R., Lartigue, G., Bénard, P., Ménard, T., Berlemont, A., 2022. A massively parallel accurate conservative level set algorithm for simulating turbulent atomization on adaptive unstructured grids. Journal of Computational Physics 458, 111075.
URL https://www.sciencedirect.com/science/article/pii/S0021999122001371

  * Janodet et al. (2019) Janodet, R., Vaudor, G., Lartigue, G., Benard, P., Moureau, V., Mercier, R., Sep. 2019. An unstructured conservative level-set algorithm coupled with dynamic mesh adaptation for the computation of liquid-gas flows. In: 29th European Conference on Liquid Atomization and Spray Systems (ILASS Europe). Paris, France.
URL https://hal.archives-ouvertes.fr/hal-02304125

  * Karypis (2003) Karypis, G., 2003. Multi-constraint mesh partitioning for contact/impact computations. In: Proceedings of the 2003 ACM/IEEE Conference on Supercomputing. SC ’03. Association for Computing Machinery, New York, NY, USA, p. 56.
URL https://doi.org/10.1145/1048935.1050206

  * Koo et al. (2014) Koo, B., Yang, J., Yeon, S., Stern, F., 2014. Reynolds and froude number effect on the flow past an interface piercing circular cylinder. Int. J. Nav. Archit. Ocean Eng. (6(3)), 529–561.
  * Kraushaar (2011) Kraushaar, M., Dec. 2011. APPLICATION OF THE COMPRESSIBLE AND LOW-MACH NUMBER APPROACHES TO LARGE-EDDY SIMULATION OF TURBULENT FLOWS IN AERO-ENGINES. Theses, Institut National Polytechnique de Toulouse - INPT.
URL https://tel.archives-ouvertes.fr/tel-00711480

  * Ling et al. (2015) Ling, Y., Zaleski, S., Scardovelli, R., 2015. Multiscale simulation of atomization with small droplets represented by a Lagrangian point-particle model. International Journal of Multiphase Flow 76, 122–143.
  * Lobovský et al. (2014) Lobovský, L., Botia-Vera, E., Castellana, F., Mas-Soler, J., Souto-Iglesias, A., 2014. Experimental investigation of dynamic pressure loads during dam break. Journal of Fluids and Structures 48, 407 – 434.
URL http://www.sciencedirect.com/science/article/pii/S0889974614000656

  * Malandain et al. (2013) Malandain, M., Maheu, N., Moureau, V., 2013. Optimization of the deflated conjugate gradient algorithm for the solving of elliptic equations on massively parallel machines. Journal of Computational Physics 238, 32–47.
URL https://www.sciencedirect.com/science/article/pii/S0021999112007280

  * Mansard and Funke (1980) Mansard, E. P., Funke, E., 1980. The measurement of incident and reflected spectra using a least squares method. In: Coastal Engineering 1980. pp. 154–172.
  * Ménard et al. (2007) Ménard, T., Tanguy, S., Berlemont, A., 2007. Coupling level set/ volume of fluid/ ghost fluid methods, validation and application to 3d simulation of the primary breakup of a liquid jet. International Journal of Multiphase Flows 33, 510–524.
  * Moureau et al. (2011) Moureau, V., Domingo, P., Vervisch, L., 2011. Design of a massively parallel cfd code for complex geometries. Comptes Rendus Mécanique 339 (2), 141 – 148, high Performance Computing.
URL http://www.sciencedirect.com/science/article/pii/S1631072110002111

  * Mukha et al. (2021) Mukha, T., Bensow, R., Liefvendahl, M., 2021. Predictive accuracy of wall-modelled large-eddy simulation on unstructured grids. Computers & Fluids 221, 104885.
URL https://www.sciencedirect.com/science/article/pii/S0045793021000517

  * Padgett et al. (2008) Padgett, J., DesRoches, R., Nielson, B., Yashinsky, M., Kwon, O.-S., Burdette, N., Tavera, E., 2008. Bridge damage and repair costs from hurricane katrina. Journal of Bridge Engineering 13 (1), 6–14.
  * Pope (1996) Pope, S., 1996. Turbulent Flows. Springer.
  * Potgieter et al. (2022) Potgieter, H., Carter, J., Henderson, D., 2022. Modeling the second harmonic in surface water waves using generalizations of nls. Water Waves 4, 23–47.
URL https://link.springer.com/article/10.1007/s42286-022-00055-7

  * Qu et al. (2021) Qu, K., Sun, W., Ren, X., Kraatz, S., Jiang, C., 2021. Numerical investigation on the hydrodynamic characteristics of coastal bridge decks under the impact of extreme waves. Journal of Coastal Research 37, 442–455.
  * Rabaud and Moisy (2013) Rabaud, M., Moisy, F., May 2013. Ship wakes: Kelvin or mach angle? Phys. Rev. Lett. 110, 214503.
URL https://link.aps.org/doi/10.1103/PhysRevLett.110.214503

  * Rapp and Melville (1990) Rapp, R., Melville, W., 1990. Laboratory measurements of deep-water breaking waves. Philosophical Transactions of the Royal Society of London 331, 735–800.
  * Santo et al. (2018) Santo, H., Taylor, P., Day, A., Nixon, E., Choo, Y., 2018. Current blockage and extreme forces on a jacket model in focussed wave groups with current. Journal of Fluids and Structures 78, 24–35.
URL https://www.sciencedirect.com/science/article/pii/S0889974617301196

  * Sutherland and O’Donoghue (1998) Sutherland, J., O’Donoghue, T., 1998. Wave phase shift at coastal structures. Journal of Waterway, Port, Coastal, and Ocean Engineering 124 (2), 90–98.
URL
https://ascelibrary.org/doi/abs/10.1061/%28ASCE%290733-950X%281998%29124%3A2%2890%29

  * Taylor and Wu (1997) Taylor, R. E., Wu, G., 1997. Interaction of steep waves with offshore structures. Philosophical Transactions: Mathematical, Physical and Engineering Sciences 355, 593–606.
  * Tryggvason et al. (2011) Tryggvason, G., Scardovelli, R., Zaleski, S., 2011. Direct Numerical Simulations of Gas-Liquid Multiphase Flows. Cambridge University Press.
  * Vaudor et al. (2017) Vaudor, G., Ménard, T., Aniszewski, W., Doring, M., Berlemont, A., 2017. A consistent mass and momentum flux computation method for two phase flows. application to atomization process. Computers & Fluids 152, 204–216.
URL https://www.sciencedirect.com/science/article/pii/S0045793017301494

  * Wernert and Favier (1999) Wernert, P., Favier, D., 1999. Considerations about the phase averaging method with application to eldv and piv measurements over pitching airfoils. Experiments in Fluids 27, 473–483.
  * Xie et al. (2017) Xie, B., Jin, P., Xiao, F., 2017. An unstructured-grid numerical model for interfacial multiphase fluids based on multi-moment finite volume formulation and THINC method. International Journal of Multiphase Flow 89, 375 – 398.
URL http://www.sciencedirect.com/science/article/pii/S0301932216302178

  * Xie et al. (2020) Xie, Z., Stoesser, T., Yan, S., Ma, Q., Lin, P., 2020. A cartesian cut-cell based multiphase flow model for large-eddy simulation of three-dimensional wave-structure interaction. Computers & Fluids 213, 104747.
URL https://www.sciencedirect.com/science/article/pii/S0045793020303170