# A Multi-Fidelity Model for Wave Energy Converters

Beatrice Battisti, Giovanni Bracco, Michel Bergmann

► To cite this version:

Beatrice Battisti, Giovanni Bracco, Michel Bergmann. A Multi-Fidelity Model for Wave Energy Converters. International Journal for Numerical Methods in Fluids, 2024, 10.1002/fld.5354. hal-04796354

HAL Id: hal-04796354

https://hal.science/hal-04796354v1

Submitted on 21 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
<figure><img src="image_1.png" /></figure>

RESEARCH ARTICLE

OPEN ACCESS

# A Multi-Fidelity Model for Wave Energy Converters

Beatrice Battisti¹,² | Giovanni Bracco¹ | Michel Bergmann²,³

¹MOREnergyLab, DIMEAS, Politecnico di Torino, Turin, Italy | ²IMB, UMR 5251, Université de Bordeaux, Talence, France | ³Inria Bordeaux Sud-Ouest, Team MEMPHIS, Talence, France

**Correspondence:** Beatrice Battisti (beatrice.battisti@inria.fr)

**Received:** 16 February 2024 | **Revised:** 12 September 2024 | **Accepted:** 26 October 2024

**Funding:** The authors acknowledge the support by the European Union's Horizon 2020 Research and Innovation Programme under the Marie Skłodowska-Curie Actions, grant agreement 872442 (ARIA).

**Keywords:** coupling methodology | multi-fidelity model | proper orthogonal decomposition | reduced order model | wave energy

## ABSTRACT

The objective of this study is to develop a three-dimensional numerical model for a floating point absorber wave energy converter in the presence of sea waves, considering its interaction with a bi-fluid flow (comprising air and water). The primary aim is to create an efficient computational tool that achieves two key objectives: firstly, reducing the computational time typically associated with high-fidelity Computational Fluid Dynamics (CFD) models, and secondly, curing the lack of accuracy of low-fidelity asymptotic or projection-based reduced-order models in regions subjected to viscous and highly nonlinear effects. To address these objectives, we propose a multi-fidelity model based on domain decomposition. This approach combines a high-fidelity CFD solver, which accurately captures the behavior in viscous and nonlinear regions, with a Reduced Order Model (ROM) based on Proper Orthogonal Decomposition (POD), tailored for weakly nonlinear regions. By integrating these components spatially, we simulate the dynamics of the floating body within a unified framework. This methodology ensures precise predictions of the body's motion for both in-sample (reproduction) and out-of-sample (prediction) configurations.

## 1 | Introduction

Numerical simulations are ubiquitous in industrial sectors such as energy, where technologies undergo a virtual cycle of concept, design, and testing before physical manufacturing. For instance, although wave energy remains an emerging sector compared to established renewable sources, it heavily relies on tractable numerical models for conception and optimization, necessitating a delicate balance between accuracy and computational costs. Wave energy converters are devices that harness kinetic and potential energy from traveling waves, converting it into mechanical or electrical energy. There are various types of converters, each based on different operational principles, deployment locations, installation types, and control systems. However, most converters generally require significant motion amplitudes to

achieve optimal efficiency. Simulating these devices is complex due to the intricate and nonlinear nature of ocean waves. This complexity is further compounded by the sophisticated control strategies, aimed at optimizing the system response.

Traditional sea-keeping numerical tools, often based on linear potential theory and commonly used in offshore structure design, are widely utilized in the initial stages of wave energy converter development, providing quick estimates of system parameters [1]. However, discrepancies persist between numerical predictions and wave tank experiments [2]. Introducing some nonlinearities to weakly or fully non-linear potential flow models partially addresses this issue [3–5], while maintaining fundamental assumptions such as inviscid flow.

<footnote>This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.</footnote>

<footnote>© 2024 The Author(s). *International Journal for Numerical Methods in Fluids* published by John Wiley & Sons Ltd.</footnote>
The utilization of high-fidelity solvers such as Computational Fluid Dynamics (CFD) or Smoothed-Particle Hydrodynamics (SPH) allows accurate simulation of all nonlinearities, but remains challenging for wave energy applications due to significant computational time requirements. Indeed, as the number of devices increases, computational costs become prohibitive, particularly in optimization analyses considering large parameter spaces and multiple queries [6, 7]. This complexity necessitates the development of accurate yet computationally efficient models.

To address this need, data-driven model reduction emerges as a strategy to extract essential features and reduce computational complexity. Reduced-Order Models (ROMs) have diverse applications spanning various fields, including thermo-acoustic problems [8], image processing [9], graphics applications [10, 11], and fluid flows [12, 13]. In the context of fluid dynamics, numerous ROMs have been developed for fluid control [14–16], interaction of flows with varying body shapes [17–19], and turbulence description [20–22]. These models find applications in various domains, ranging from fluid-structure interaction in cardiovascular problems [23] to the numerical simulation of the wake behind a wind turbine [24], and the depiction of the dynamics of multi-phase flows [25], as in the case of this study. However, the application of ROMs to highly nonlinear problems like wave energy converters necessitates careful consideration.

Among the model reduction techniques, the Proper Orthogonal Decomposition (POD) method emerges as an effective means for deriving a reduced basis in high-dimensional flow systems. Intrusive ROMs, typically derived through POD and Galerkin projection methods or (least square) Petrov-Galerkin approaches [26, 27], are attractive for linear and weakly nonlinear model reductions and find success in various research domains. For instance, some studies involve different phases such as fluid and solid [28, 29], heat transfer [30], and porous media [31]. However, Galerkin ROMs might encounter challenges related to instability, and the use of the linear affine POD representation can limit efficiency for highly nonlinear problems [32], even for Petrov-Galerkin projections. Quadratic POD approximations [33] and fully nonlinear approximations [34] have been recently introduced to mitigate the Kolmogorov barrier in nonlinear projection ROMs. However, their intrusive nature often necessitates complex coding modifications [35]. The complexity of intrusive projection-based ROMs is particularly evident when the computational source code is unavailable, as is often the case with some commercial software. To overcome this limitation, we propose a non-intrusive ROM approach where POD ROM is obtained in a Galerkin-free way, coupled with a high-fidelity solver implemented in commercial code.

Similar to other coupling methodologies [36, 37], this approach relies on domain decomposition, for which ROMs are particularly effective [38–40]. In this study, the computational domain is decomposed into low-fidelity (POD) and high-fidelity (Navier-Stokes solver) subdomains. The temporal coefficients of the POD expansion are obtained through a minimization problem defined on both an overlapping domain between the two subdomains and on a domain (discrete points defining sensors) where incoming flow is actually imposed. The novelty of this paper lies in applying the established

multi-fidelity, component-based approach based on POD, to a three-dimensional, multi-phase flow with an evolving free surface. Additionally, the use of sensors outside the high-fidelity domain to propagate information elevates the POD model from a conventionally passive role. To the best of the authors' knowledge, these techniques have not been documented in the literature yet. Specifically, the proposed method may have diverse potential applications, and for this study, it has been applied to the case of wave energy, for which the high-fidelity domain is chosen to be a small surrounding of the floating body (the wave energy converter), where highly nonlinear phenomena usually occur, and the other POD domain is used to both give boundary conditions to the high-fidelity solver, and to propagate the information elsewhere using global POD modes. This multi-fidelity numerical model is used to predict the temporal evolutions of velocity field, the water surface, and the body dynamics in wave energy converter simulations.

The paper is organized as follows: The multi-fidelity model, which couples the high-fidelity Navier-Stokes model and the low-fidelity POD model, is presented in Section 2. A in-depth sensitivity analysis of several parameters involved in the multi-fidelity model, such as the number of POD modes or the size of the different subdomains, is described for a simple flow configuration in Section 3. A floating point absorber wave energy converter is then studied in Section 4. Finally, some conclusions are drawn in Section 5.

## 2 | Numerical Multi-Fidelity Model

### 2.1 | Flow Configuration

The objective is to characterize the performance of a wave energy harvesting device under the influence of sea waves. A three-dimensional domain is considered, and a two-dimensional vertical section is schematically represented in Figure 1. We consider a point absorber-type structure, wherein the sole degree of freedom is the vertical motion. The densities and viscosities are represented by $\rho$ and $\mu$, with subscripts *a* and *w* indicating air and water, respectively. The floater, characterized by a density $\rho_s$,

<figure><img src="image_2.png" /><figcaption>**FIGURE 1** | Sketch of the flow configuration. [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>
occupies domain $\Omega_s$, while the fluid (comprising air and water) is flowing within domain $\Omega_f$. The fluid domain $\Omega_f$ is further partitioned into two distinct regions, $\Omega_a$ for air and $\Omega_{iw}$ for water, with an interface denoted by $\Gamma_f$.

Given the nature of a floating structure, we have $\rho_a < \rho_s < \rho_{iw}$.

Before introducing the multi-fidelity model (Section 2.4), we give a description of the two models involved in the coupling procedure, the high-fidelity model (Section 2.2) and the low-fidelity POD approximation (Section 2.3).

## 2.2 | High-Fidelity Numerical Model

The flow described in the previous section is modelled by the two-phase Navier-Stokes equations. Let $\alpha$ represent the local fraction of water, referred to as the Volume Of Fluid (VOF), where $0 \le \alpha \le 1$ [41]. The fluid density and dynamic viscosity can be expressed as follows:

$$ \rho(\mathbf{x}, t) = \rho_a + (\rho_{iw} - \rho_a) \alpha(\mathbf{x}, t) \quad (1) $$

$$ \mu(\mathbf{x}, t) = \mu_a + (\mu_{iw} - \mu_a) \alpha(\mathbf{x}, t) \quad (2) $$

where $\alpha = 0$ and $\alpha = 1$ denote air and water, respectively. The bi-fluid interface is approximated by $\alpha = 0.5$.

Let $D(\mathbf{u}) = \frac{\nabla u + (\nabla u)^T}{2}$ denote the viscous part of the stress tensor, where the velocity field is $\mathbf{u} = (u, v, w)$ with $u, v$, and $w$ representing the velocity components in the $x, y$, and $z$ directions, respectively. Let also $p$ be the pressure field and $g$ the gravity acceleration along the $y$ axis. The incompressible Navier-Stokes equations for both fluids in the domain $\Omega_f = \Omega_a \cup \Omega_{iw}$ are:

$$ \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla)\mathbf{u} = -\frac{1}{\rho}\nabla p + \frac{1}{\rho}\nabla \cdot 2\mu D(\mathbf{u}) + g \quad \text{in } \Omega_f \quad (3) $$

$$ \nabla \cdot \mathbf{u} = 0 \quad \text{in } \Omega_f \quad (4) $$

where the temporal evolution of the bi-fluid interface, and thus $\mu$ and $\rho$, is governed by the transport equation:

$$ \frac{\partial \alpha}{\partial t} + \mathbf{u} \cdot \nabla \alpha = 0 \quad \text{in } \Omega_f \quad (5) $$

This system is complemented with initial conditions $\mathbf{u}(\mathbf{x}, t=0) = \mathbf{u}_0$, $p(\mathbf{x}, t=0) = p_0$, and $\alpha(\mathbf{x}, t=0) = \alpha_0$, boundary conditions on the external boundary $\Gamma_e$ for velocity $\mathbf{u}(\mathbf{x}, t)$ and VOF $\alpha(\mathbf{x}, t)$, and boundary conditions on the structure boundary $\Gamma_s$, as well as conditions through the bi-fluid interface $\Gamma_f$ [42], defined at $\alpha = 0.5$, and whose evolution follows Equation (5).

The motion of the floating rigid body is governed by Newton's laws:

$$ m \frac{d\bar{\mathbf{u}}}{dt} = F_{\text{ext}} \quad (6) $$

$$ \frac{dJ\omega}{dt} = M_{\text{ext}} \quad (7) $$

where $m$ and $J$ are the mass and inertia matrix of the body, and $\bar{\mathbf{u}}$ and $\omega$ denote the linear and angular velocities.

The forces and torques are computed by:

$$ F_{\text{ext}} = - \int_{\Gamma_s} \mathbb{T}(\mathbf{u}, p) \, \mathbf{n} \, d\mathbf{x} + g \quad (8) $$

$$ M_{\text{ext}} = - \int_{\Gamma_s} \mathbf{r} \wedge \mathbb{T}(\mathbf{u}, p) \, \mathbf{n} \, d\mathbf{x} \quad (9) $$

where $\mathbb{T}(\mathbf{u}, p) = -p\mathbf{I} + \mu(\nabla\mathbf{u} + \nabla\mathbf{u}^T)$ is the stress tensor, $\mathbf{n}$ is the unit outward vector to $\Gamma_s$, and $\mathbf{r} = \mathbf{x} - \mathbf{x}_G$ with $\mathbf{x}_G$ representing the body center of mass.

In this study, the numerical solution of the Navier-Stokes equations is computed using the commercial software StarCCM+ [43]. The governing equations are discretized on an overset Chimera mesh, where a foreground body-fitted grid refined near the body is superimposed on an octree backward grid refined around the air-water interface. A typical computational mesh used in this study is shown in Figure 2.

For open flows, artificial boundary conditions are typically applied at external boundaries. These may include Dirichlet or Neumann conditions, or more sophisticated absorbing or non-reflecting conditions [44, 45], depending on the complexity of the problem under consideration. Imposing appropriate boundary conditions remains one of the main challenges in Computational Fluid Dynamics (CFD). In open flows, such as the one described in Figure 1, boundary conditions are usually applied far from the body, resulting in a large computational domain. The computational costs required to solve the Navier-Stokes equations in such a computational domain can become prohibitively large. To decrease these costs, reduced-order models can be employed, either as a substitution model or to provide more suitable boundary conditions that allow for a reduction in the size of the computational domain [46, 47].

Instead of computing the solution of the Navier-Stokes equations on a fine mesh with possibly millions of cells, we can approximate the solution using a few global basis functions (typically on the order of $O(10)$). In this study, the global basis functions are computed using Proper Orthogonal Decomposition (POD).

<figure><img src="image_3.png" /><figcaption>FIGURE 2 | Sketch of the computational mesh. [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>
## 2.3 | Low-Fidelity Subspace: Proper Orthogonal Decomposition

The goal is to approximate the numerical high-fidelity solution $U_h$ of the governing Equations (3–5) by a POD ansatz $\bar{U}$ with few $N_r = \mathcal{O}(10)$ global basis functions $\{\phi_i\}_{i=1}^{N_r}$, that is

$$U_h(\mathbf{x}, t) \approx \tilde{U}(\mathbf{x}, t) = \bar{U}(\mathbf{x}) + \sum_{n=1}^{N_r} a_n(t)\phi_n(\mathbf{x}) \quad (10)$$

The field $U$ can be the velocity field $u$, the pressure field $p$, the density field $\rho$, the VOF field $\alpha$, or even a mixed field $(u, p, \rho, \alpha)$ (or other alternatives) [48], and $\bar{U}$ is a reference field. The latter is usually the average in time, a shift or non equilibrium mode [49, 50], or zero. When the field $U$ is combined, the different variables may have varying orders of magnitude. In such instances, weights can be assigned to the variables to balance their significance or to emphasize their importance relative to the other components. In this paper, no weights have been applied, or equivalently, all weights are set to unity.

The Proper Orthogonal Decomposition (POD) was first introduced in turbulence by Lumley [51] in 1967 to identify the coherent structures in a turbulent flow. For a review of the POD, readers are referred to [52, 53]. The goal is to find the deterministic functions $\phi(\mathbf{x})$ that are most similar in an average sense to data $U(\mathbf{x}, t) = U_h(\mathbf{x}, t) - \bar{U}(\mathbf{x})$. In what follows, the snapshot method introduced by Sirovich [54] is used.

The construction of the low-fidelity subspace, the POD basis functions, is called the offline stage. It is based on the following steps.

First, one has to collect a set of $N_s$ snapshots $\{U_h(\mathbf{x}, t_j)\}_{j=1}^{N_s}$ where $t_j = T_0 + (j-1)\Delta T$, for given values of $T_0$ and $\Delta T$. Note that snapshots can also be sampled on different input parameters $\eta_i$ (as it will be the case hereafter for different wave parameters). Without loss of generality, the POD description only considers time, and we then build $U(\mathbf{x}, t_j) = U_h(\mathbf{x}, t_j) - \bar{U}(\mathbf{x})$, where $\bar{U}(\mathbf{x})$ is the temporal average of the $N_s$ snapshots $U_h(\mathbf{x}, t_j)$.

In what follows, the snapshots are extracted (and possibly interpolated) on a uniform Cartesian grid composed by $n_c$ cells. In this study, we consider $U = (u, v, w, \alpha)^T$. Separating the POD modes for different variables would allow for the definition of a specific number of modes for each variable. Conversely, a potential drawback of computing POD basis functions from velocity and VOF fields together is that it requires the same number of modes for both, which may result in sub-optimal data compression. However, the overall results from the coupled simulation using separated POD modes were either worse or comparable to those using global modes, with stability issues arising. It seems that the necessary interpolation during the information exchange

mitigated any improvements, suggesting that a strong connection between the velocity and the VOF is required to stabilize the boundary conditions. Additionally, the main advantage of the global approach is that only one POD is needed, ensuring all variables remain closely linked. If needed, one could also add the (scaled) pressure field to $U$.

The snapshots matrix is $S = [U(\mathbf{x}, t_1)^T \ U(\mathbf{x}, t_2)^T \ \dots \ U(\mathbf{x}, t_{N_s})^T] \in \mathbb{R}^{4 \times N_s}$. Numerically, we have:

$$S = \begin{pmatrix} u(x_1, t_1) & u(x_1, t_2) & \cdots & u(x_1, t_{N_s}) \\ \vdots & \vdots & & \vdots \\ u(x_{n_c}, t_1) & u(x_{n_c}, t_2) & \cdots & u(x_{n_c}, t_{N_s}) \\ v(x_1, t_1) & v(x_1, t_2) & \cdots & v(x_1, t_{N_s}) \\ \vdots & \vdots & & \vdots \\ v(x_{n_c}, t_1) & v(x_{n_c}, t_2) & \cdots & v(x_{n_c}, t_{N_s}) \\ w(x_1, t_1) & w(x_1, t_2) & \cdots & w(x_1, t_{N_s}) \\ \vdots & \vdots & & \vdots \\ w(x_{n_c}, t_1) & w(x_{n_c}, t_2) & \cdots & w(x_{n_c}, t_{N_s}) \\ \alpha(x_1, t_1) & \alpha(x_1, t_2) & \cdots & \alpha(x_1, t_{N_s}) \\ \vdots & \vdots & & \vdots \\ \alpha(x_{n_c}, t_1) & \alpha(x_{n_c}, t_2) & \cdots & \alpha(x_{n_c}, t_{N_s}) \end{pmatrix} \in \mathbb{R}^{4 \times N_s} \quad (11)$$

The temporal correlation matrix is

$$R = S^T S \in \mathbb{R}^{N_s \times N_s} \quad (12)$$

If one wants to consider another parameter $\eta$, the snapshot matrix becomes $S = (S_{\eta_1}, S_{\eta_2}, \dots, S_{\eta_{N_m}})$, where $N_m$ is the number of sampling points for parameters $\eta$, and $S_{\eta_i}$ is the snapshot matrix (11) built with $\eta = \eta_i$, and we finally have $N_s := N_s \times N_m$ snapshots.

We then solve the following eigenvalues problem

$$R\psi_n = \lambda_n\psi_n, \quad n = 1, \dots, N_s \quad (13)$$

where the $\lambda_1 \ge \lambda_2 \ge \dots \ge \lambda_{N_s} \ge 0$.

The POD basis functions, also called POD modes, are:

$$\phi_n = \frac{1}{\sqrt{\lambda_n}} S \psi_n, \quad n = 1, \dots, N_s \quad (14)$$

and the associated POD temporal coefficients are $a_n = \frac{\psi_n}{\sqrt{\lambda_n}}$. From this point forward, these temporal coefficients will also be referred to as projection coefficients, $a^{\text{proj}}$.

Finally, from the relative information content $\text{RIC}(n) = \sum_{i=1}^n \lambda_i / \sum_{i=1}^{N_s} \lambda_i$, the POD basis is truncated to keep only the $N_r \ll N_s$ most dominant modes $\{\phi_i\}_{i=1}^{N_r}$, such that $\text{RIC}(N_r) > \delta$, with a user defined threshold $\delta$ usually taken to be larger than 99%.

Let $\Phi \in \mathbb{R}^{4 \times N_r}$ be the matrix built with the $N_r$ first POD modes, that is
$$ \mathbf{\Phi} = \begin{pmatrix} \phi_1^u(x_1) & \phi_2^u(x_1) & \cdots & \phi_{N_r}^u(x_1) \\ \vdots & \vdots & & \vdots \\ \phi_1^u(x_{n_c}) & \phi_2^u(x_{n_c}) & \cdots & \phi_{N_r}^u(x_{n_c}) \\ \phi_1^v(x_1) & \phi_2^v(x_1) & \cdots & \phi_{N_r}^v(x_1) \\ \vdots & \vdots & & \vdots \\ \phi_1^v(x_{n_c}) & \phi_2^v(x_{n_c}) & \cdots & \phi_{N_r}^v(x_{n_c}) \\ \phi_1^w(x_1) & \phi_2^w(x_1) & \cdots & \phi_{N_r}^w(x_1) \\ \vdots & \vdots & & \vdots \\ \phi_1^w(x_{n_c}) & \phi_2^w(x_{n_c}) & \cdots & \phi_{N_r}^w(x_{n_c}) \\ \phi_1^a(x_1) & \phi_2^a(x_1) & \cdots & \phi_{N_r}^a(x_1) \\ \vdots & \vdots & & \vdots \\ \phi_1^a(x_{n_c}) & \phi_2^a(x_{n_c}) & \cdots & \phi_{N_r}^a(x_{n_c}) \end{pmatrix} \in \mathbb{R}^{4 \times N_r} \quad (15) $$

The Equation (10) writes:

$$ U \approx \tilde{U} = \bar{U} + \Phi a \quad (16) $$

with $a = (a_1, a_2, \dots, a_{N_r})^T$.

Defining the sub-matrices

$$ \begin{align} \Phi^u &= \begin{pmatrix} \phi_1^u(x_1) & \phi_2^u(x_1) & \cdots & \phi_{N_r}^u(x_1) \\ \vdots & \vdots & & \vdots \\ \phi_1^u(x_{n_c}) & \phi_2^u(x_{n_c}) & \cdots & \phi_{N_r}^u(x_{n_c}) \end{pmatrix} \in \mathbb{R}^{n_c \times N_r} \notag \\ \Phi^v &= \begin{pmatrix} \phi_1^v(x_1) & \phi_2^v(x_1) & \cdots & \phi_{N_r}^v(x_1) \\ \vdots & \vdots & & \vdots \\ \phi_1^v(x_{n_c}) & \phi_2^v(x_{n_c}) & \cdots & \phi_{N_r}^v(x_{n_c}) \end{pmatrix} \in \mathbb{R}^{n_c \times N_r} \notag \\ \Phi^w &= \begin{pmatrix} \phi_1^w(x_1) & \phi_2^w(x_1) & \cdots & \phi_{N_r}^w(x_1) \\ \vdots & \vdots & & \vdots \\ \phi_1^w(x_{n_c}) & \phi_2^w(x_{n_c}) & \cdots & \phi_{N_r}^w(x_{n_c}) \end{pmatrix} \in \mathbb{R}^{n_c \times N_r} \notag \\ \Phi^\alpha &= \begin{pmatrix} \phi_1^\alpha(x_1) & \phi_2^\alpha(x_1) & \cdots & \phi_{N_r}^\alpha(x_1) \\ \vdots & \vdots & & \vdots \\ \phi_1^\alpha(x_{n_c}) & \phi_2^\alpha(x_{n_c}) & \cdots & \phi_{N_r}^\alpha(x_{n_c}) \end{pmatrix} \in \mathbb{R}^{n_c \times N_r} \tag{17} \end{align} $$

each physical field can be approximated onto these bases:

$$ u \approx \tilde{u} = \bar{u} + \Phi^u a $$

$$ v \approx \tilde{v} = \bar{v} + \Phi^v a $$

$$ w \approx \tilde{w} = \bar{w} + \Phi^w a $$

$$ \alpha \approx \tilde{\alpha} = \bar{\alpha} + \Phi^\alpha a $$

As previously mentioned, one advantage of computing a global POD basis with $U = (u, v, w, \alpha)^T$ is that the same temporal coefficients $a$ are used for each component.

Note that for out-of-sample computation, that is for $t \neq t_j$ or $\eta \neq \eta_i$ (different operating conditions), the temporal coefficients $a$ have to be computed using an external model. Usually, a projection-based reduced order model is used, like a Galerkin or Petrov-Galerkin model. However, for the bi-fluid case studied here, the last term in Equation (3) would lead to a fourth-order

<figure><img src="image_4.png" /><figcaption>FIGURE 3 | Domain decomposition for the coupling methodology. The overlapping domain $\Omega_o$ is dark-gray colored. The high-fidelity domain $\Omega_{hf}$ is the union of the white and dark gray domains, and the low-fidelity domain $\Omega_{lf}$ is the union of the light and dark gray domains.</figcaption></figure>

tensor that is not efficient for model order reduction [55]. Hyper-reduction techniques [56, 57] can be in principle used to mitigate that problem. However, these intrusive approaches are impossible in the context of commercial software where the sources, and thus the numerical residuals, are not available. A modified and extended version of the Galerkin-free approach introduced in [46] is used and presented hereafter.

In what follows, the POD basis functions are numerically computed on a Cartesian mesh that corresponds to the coarsest cell size in the high-fidelity mesh (see Figure 2). If a body is present in the computational domain, mesh points located within all potential positions of the body during its motion are removed—typically by excluding a Cartesian box. One approach is to compute the POD basis functions only on the union of the two gray domains in Figure 3, as the body moves exclusively within the white domain.

## 2.4 | Multi-Fidelity Hybrid Model

### 2.4.1 | General Process

The multi-fidelity model proposed here couples a high-fidelity solver with a low-fidelity solver, where the Proper Orthogonal Decomposition (POD) temporal coefficients ${a_i}^{N_r}_{i=1}$ are obtained using localized information. This information is provided by high-fidelity solutions on subdomains, and possibly by imposed values at external local sensors. Let $\Omega_{lf}$ denote the subdomain where the low-fidelity model is used, $\Omega_{hf}$ denote the subdomain where the high-fidelity model is used, and $\Omega_o = \Omega_{lf} \cap \Omega_{hf}$ represent the overlapping domain where localized information is exchanged between high-fidelity and low-fidelity models. Let ${x_i^S}^M_{i=1}$ be the locations of $M$ sensors, where some information $\hat{U}_i = (\hat{u}_i, \hat{v}_i, \hat{w}_i, \hat{\alpha}_i)^T$ is given (at least partially). For some applications, we may have $M=0$. A sketch of the domain decomposition is provided in Figure 3.
The main requirements for the high-fidelity domain are as follows:

* It should be as small as possible to reduce computational costs.

* It should be located in regions where the POD representation is not accurate enough, such as highly nonlinear regions, especially nearby structures like floating bodies.

* It should surround all possible unsteady locations of the floater (if any).

* It should have at least $N_r - M$ discrete points to obtain a well-posed minimization problem for ${a_i}_{i=1}^{N_r}$ described hereafter.

The objective is to find the temporal coefficients $a^*(t_j) = (a_1^*, a_2^*, \dots, a_{N_r}^*)^T$, at a given time $t_j$, where $t_j$ represents the element of a set of discrete time, such that the reduced solution $\tilde{U}(\mathbf{x}, t_j) = \overline{U}(\mathbf{x}) + \Phi(\mathbf{x}) a(t_j)$ is as close as possible to the high-fidelity solution $U_{hf}(\mathbf{x}, t_j)$ on the overlapping domain $\Omega_o$, while also being close to the observed (and possibly imposed) values ${\hat{U}_i}_{i=1}^M(\mathbf{x}, t_j)$ at the $M$ sensor points. Similarly to data $\tilde{U}$, the field $\tilde{U}$ evaluated at a point-wise location $\mathbf{x}_i^S$ is denoted $\tilde{U}_i$. Depending on the sensor locations, interpolations may be required to compute $\tilde{U}_i$.

Mathematically, at a given time $t_j$, the temporal coefficients $a^*(t_j)$ satisfy:

$$ a^* = \arg \min_{a} \left( \frac{1}{|\Omega_o|} \int_{\Omega_o} \left( U_{hf} - \tilde{U} \right)^2 d\mathbf{x} + \frac{\beta}{M} \sum_{i=1}^{M} \left( \hat{U}_i - \tilde{U}_i \right)^2 \right) \\ \text{with } \tilde{U} = \overline{U} + \Phi a \quad (18) $$

where the parameter $\beta$ is used to weight the sensor contributions, and the space-time dependencies have been omitted for the sake of readability. For low values of $\beta$, the Reduced Order Model (ROM) approximation, $\tilde{U} = \overline{U} + \Phi a^*$, is mostly influenced by the Full Order Model (FOM) solution, $U_{hf}$. For large values of $\beta$, the ROM is mostly influenced by the sensor data, $\tilde{U}$, thus providing new information to the FOM. Values in between, $\beta \approx 1$, offer a compromise between FOM and sensor data contributions, which are well balanced in Equation (18). The value of $\beta$ largely depends on the specific problem and should be assessed accordingly.

In what follows, we consider only the Volume of Fluid (VOF), $\hat{a}_i$, as measurable by sensors. These data represent water elevation, which is usually easy to measure. The minimization problem (18) is thus:

$$ a^* = \arg \min_a \left( \frac{1}{|\Omega_o|} \int_{\Omega_o} \left( U_{hf} - \tilde{U} \right)^2 d\mathbf{x} + \frac{\beta}{M} \sum_{i=1}^{M} \left( \hat{a}_i - \tilde{a}_i \right)^2 \right) \quad (19) $$

However, the following minimization problem:

$$ a^* = \arg \min_a \left( \frac{1}{|\Omega_o|} \int_{\Omega_o} \left( \alpha_{hf} - \tilde{\alpha} \right)^2 d\mathbf{x} + \frac{\beta}{M} \sum_{i=1}^{M} \left( \hat{a}_i - \tilde{a}_i \right)^2 \right) \quad (20) $$

with $\tilde{\alpha} = \bar{\alpha} + \Phi^a a$, provides better results and is used in the remaining part of this study. This choice is driven not only by the

aim for consistency in the minimization equation, but also by its proven ability to yield the most stable and accurate results. This is likely because it avoids conflicting minimization outcomes for the different variables, and mitigates inevitable interpolation filtering effects. For further details, the interested reader may refer to the doctoral thesis [55].

Each physical field used in the POD snapshots can finally be approximated as:

$$ \tilde{u} = \bar{u} + \Phi^u a^* \quad (21) $$

$$ \tilde{v} = \bar{v} + \Phi^v a^* \quad (22) $$

$$ \tilde{w} = \bar{w} + \Phi^w a^* \quad (23) $$

$$ \tilde{\alpha} = \bar{\alpha} + \Phi^\alpha a^* \quad (24) $$

## 2.4.2 | Continuous Multi-Fidelity Model

The multi-fidelity model finally writes:

$$ \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla)\mathbf{u} = -\frac{1}{\rho}\nabla p + \frac{1}{\rho}\nabla \cdot 2\mu D(\mathbf{u}) + g \quad \text{in } \Omega_{hf} \quad (25) $$

$$ \nabla \cdot \mathbf{u} = 0 \quad \text{in } \Omega_{hf} \quad (26) $$

$$ \frac{\partial \alpha}{\partial t} + \nabla \cdot (\alpha \mathbf{u}) = 0 \quad \text{in } \Omega_{hf} \quad (27) $$

$$ u = \tilde{u} = \bar{u} + \Phi^u a^* \quad \text{on } \Gamma_{hf} \quad (28) $$

$$ v = \tilde{v} = \bar{v} + \Phi^v a^* \quad \text{on } \Gamma_{hf} \quad (29) $$

$$ w = \tilde{w} = \bar{w} + \Phi^w a^* \quad \text{on } \Gamma_{hf} \quad (30) $$

$$ \alpha = \tilde{\alpha} = \bar{\alpha} + \Phi^\alpha a^* \quad \text{on } \Gamma_{hf} \quad (31) $$

$$ a^* = \arg \min_a \left( \frac{1}{|\Omega_o|} \int_{\Omega_o} \left( \alpha_{hf} - \tilde{\alpha} \right)^2 d\mathbf{x} + \frac{\beta}{M} \sum_{i=1}^{M} \left( \hat{a}_i - \tilde{a}_i \right)^2 \right) \\ \text{with } \tilde{\alpha} = \bar{\alpha} + \Phi^a a \quad (32) $$

For the sake of readability, the space and time dependencies are omitted, but it should be noted that all variables (velocity, density, pressure, VOF) depend on both space and time. The POD modes are defined in space, while the temporal coefficients vary over time. The system of Equations (25–32) is discretized in space across the different domains using meshes, and is solved at specific, discrete time instants.

Here, Dirichlet boundary conditions are used for both velocity and the VOF fields. Although there may be various options for the boundary conditions, imposing Dirichlet boundary conditions on all boundaries of the domain provides the best solution, given the constraints of the commercial software used.

Note that the field reconstructions (28-31) actually extend over the entire domain $\Omega_{hf}$ (including the boundaries $\Gamma_{hf}$), thus providing a propagation model.
### 2.4.3 | Discrete Approach

The high-fidelity and low-fidelity solutions can be computed on different meshes. Generally, the high-fidelity simulations can be computed on any type of mesh. Concerning the low-fidelity solution, the only requirement is a constant, time-invariant mesh, since the POD modes are tied to a specific spatial distribution. Interpolations are, thus, necessary to bridge the differences between the potentially very different computational frameworks. In this paper, Cartesian meshes are used for both CFD and POD computations. The high-fidelity solution $U$ is computed on an octree grid refined near the interface, possibly coupled with an overset grid (when floating structures are considered). The low-fidelity solution $\tilde{U}$ is computed on a uniform Cartesian grid, corresponding to the coarsest level of refinement in the octree grid. Interpolations are thus required for Equations (28–32), where information between subdomains has to be exchanged. Different interpolation functions are introduced:

* $I_↓$ is an interpolator from the high-fidelity mesh (octree) to the low-fidelity mesh (Cartesian uniform). It is used, for instance, to numerically compute the integral over $\Omega_o$ (the first term in (32)).

* $I_↑$ is an interpolator from the low-fidelity mesh to the high-fidelity mesh. It is used, for instance, to evaluate the Dirichlet boundary conditions (28–31) on $\Gamma_{hf}$, $u = I_↑(\tilde{u})$ (similarly for $v, w$, and $\alpha$). The property $I_↑ = I_↓^{-1}$ is expected, however, it might not be satisfied numerically.

* $I_s$ is an interpolator from the low-fidelity Cartesian mesh to sensor locations $\mathbf{x}_i^S \in \Omega_{lf}$, that is $\tilde{\alpha}_i(\mathbf{x}_i^S) = I_s(\tilde{\alpha}_i(\{\mathbf{x}_i^{lf}\}))$, where $\{\mathbf{x}_i^{lf}\}$ denotes formally the low fidelity Cartesian mesh stencil used to interpolate on point location $\mathbf{x}_i^S$. It is used to evaluate the second term in (32).

Equations (28–31) become:

$$u = I_↑(\tilde{u}) \quad v = I_↑(\tilde{v}) \quad w = I_↑(\tilde{w}) \quad \alpha = I_↑(\tilde{\alpha}) \quad \text{on } \Gamma_{hf} \quad (33)$$

and Equation (32) becomes:

$$\mathbf{a}^* = \arg\min_{\mathbf{a}} \left( \frac{1}{|\Omega_o|} \int_{\Omega_o} (I_↓(\alpha_{hf}) - \tilde{\alpha})^2 \mathrm{d}\mathbf{x} + \frac{\beta}{M} \sum_{i=1}^{M} (\hat{\alpha}_i - I_s(\tilde{\alpha}_i))^2 \right) \quad (34)$$

If, alternatively, the high-fidelity octree mesh was used for the numerical computation of the integral over the overlapping domain $\Omega_o$, the first term in Equation (34) would be replaced by $\frac{1}{|\Omega_o|} \int_{\Omega_o} (\alpha_{hf} - I_↑(\tilde{\alpha}))^2 \mathrm{d}\mathbf{x}$. However, this paper employs the low-fidelity Cartesian mesh, as it remains static over time, ensuring a one-to-one correspondence for all mesh points.

### 2.4.4 | Errors Computations

In what follows, $\xi$ can be a velocity component $u, v$ or $w$, or the VOF $\alpha$. Let $\xi_{HF}$ be the reference high-fidelity solution computed on the large computational domain $\Omega_{HF}$. The two following fields are defined from the coupling algorithm: the field $\xi_{hf}$ is the high-fidelity solution computed on the small domain $\Omega_{hf}$, and the

field $\tilde{\xi}$ is the low-fidelity solution computed in domain $\Omega_{lf}$. The relative errors are:

$$\epsilon_{hf}^{\tilde{\xi}}(t) = \frac{\sqrt{\int_{\Omega_{hf}} (\xi_{HF}(\mathbf{x}, t) - \xi_{hf}(\mathbf{x}, t))^2 \mathrm{d}\mathbf{x}}}{\sqrt{\int_{\Omega_{hf}} \xi_{HF}(\mathbf{x}, t)^2 \mathrm{d}\mathbf{x}}} \quad (35)$$

$$\epsilon_{lf}^{\tilde{\xi}}(t) = \frac{\sqrt{\int_{\Omega_{lf} \setminus \Omega_o} (\xi_{HF}(\mathbf{x}, t) - \tilde{\xi}(\mathbf{x}, t))^2 \mathrm{d}\mathbf{x}}}{\sqrt{\int_{\Omega_{lf} \setminus \Omega_o} \xi_{HF}(\mathbf{x}, t)^2 \mathrm{d}\mathbf{x}}} \quad (36)$$

Solutions $\xi_{HF}$, $\xi_{hf}$ and $\tilde{\xi}$ are computed on different grids. It is, thus, more convenient to compute the error on the same grid, and we chose the uniform Cartesian grid used to computed the low-fidelity solution. Numerically, errors (35) and (36) are approximated by:

$$\tilde{\epsilon}_{hf}^{\tilde{\xi}}(t) = \frac{\| I_↓(\xi_{HF}(t)) - I_↓(\xi_{hf}(t)) \|_{\Omega_{hf}}}{\| I_↓(\xi_{HF}(t)) \|_{\Omega_{hf}}} \quad (37)$$

$$\tilde{\epsilon}_{lf}^{\tilde{\xi}}(t) = \frac{\| I_↓(\xi_{HF}(t)) - \tilde{\xi}(t) \|_{\Omega_{lf} \setminus \Omega_o}}{\| I_↓(\xi_{HF}(t)) \|_{\Omega_{lf} \setminus \Omega_o}} \quad (38)$$

where $\| \cdot \|_{\Omega_i}$ stands for the discrete $L_2$-norm on a domain $\Omega_i$.

The total numerical error on the whole domain $\Omega_f$ can be approximated by $\tilde{\epsilon}_{hf}^{\tilde{\xi}}(t) = \tilde{\epsilon}_{hf}^{\tilde{\xi}}(t) + \tilde{\epsilon}_{lf}^{\tilde{\xi}}(t)$. However, note that $\tilde{\epsilon}_{hf}^{\tilde{\xi}}(t)$ is not equivalent to the relative $L_2$ error on $\Omega_{HF}$. In what follows, we will only consider temporal average error over a given horizon $T$, and we define:

$$\tilde{\epsilon}_{hf}^{\tilde{\xi}} = \int_{T} \tilde{\epsilon}_{hf}^{\tilde{\xi}}(t) \mathrm{d}t \quad (39)$$

$$\tilde{\epsilon}_{lf}^{\tilde{\xi}} = \int_{T} \tilde{\epsilon}_{lf}^{\tilde{\xi}}(t) \mathrm{d}t \quad (40)$$

## 3 | Parameter Sensitivity Analysis

As shown in the previous sections, various parameters must be selected, such as the number of Proper Orthogonal Decomposition (POD) modes or the sizes of the computational domains, $\Omega_{hf}$ and $\Omega_o$.

In this analysis, we consider a simple dam break-like configuration without a body nor wave forcing ($M = 0$). Introducing a body or an incoming wave would affect the overall simulation, and slightly influence the sensitivity analysis. The test case chosen is computationally affordable and allows for in-depth parameter sensitivity analysis. While other flow configurations may yield different optimal parameter values, this study serves to identify the parameters crucial to the coupling strategy.

A three-dimensional domain $\Omega_{HF}$ with dimensions $[L_X \times L_Y \times L_Z] = [10 \times 6 \times 0.4]$ m is considered. The small size of $L_Z$ ensures reasonable computational costs.

The initial position of the water-air interface $\Gamma_f$ is analytically defined by:

$$\Gamma_f^0 = H_w e^{-\frac{1}{2}(x+\frac{3}{2}y)^2} + H_s \quad (41)$$
where parameter $H_\omega$ determines the initial wave height relative to the sea level $H_s$ at $x \to \pm\infty$. The height $H_s$ is adjusted to achieve a desired water volume in the computational domain. Here, we consider $H_s = 0$ m, and $0 \text{ m} \le H_\omega \le 1$ m to obtain smooth waves commonly observed in the wave energy sector, where extreme waves are rare.

The overall configuration is depicted in Figure 4.

<figure><img src="image_5.png" /><figcaption>**FIGURE 4** | Sketch of the flow configuration with different computational domains. The initial condition (41) for $H_\omega = 0.6$ m and $H_s = 0$ m is shown by the blue line. [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

No-slip boundary conditions, $\mathbf{u} = \mathbf{0}$ m/s, are enforced on the boundaries orthogonal to the $x$ and $y$ directions, while symmetry boundary conditions are applied to the boundary orthogonal to the $z$ direction. A few snapshots of the CFD solution across a two-dimensional section of the domain $\Omega_{HF}$ are shown in Figure 5.

The objective of this section is to conduct a sensitivity analysis of various physical and numerical parameters. Numerical sensitivities pertain to mesh discretization for the Full Order Model (FOM) (Section 3.1), temporal coupling scheme (Section 3.2), number of snapshots $N_s$, and number of POD modes $N_r$ (Section 3.3). Geometrical sensitivities involve the size of the small high-fidelity domain $\Omega_{hf}$ (Section 3.4) and the size of the overlapping domain $\Omega_o$ (Section 3.5).

Sensitivity analyses are conducted independently on each parameter while keeping all others fixed, potentially introducing biases in the analysis. To mitigate such biases, ideal scenarios will be considered for fixing the remaining parameters.

## 3.1 | Mesh Sensitivity Analysis

The primary objective is to compute an accurate solution $U_h$ of the Navier-Stokes Equations (3-5) over the entire domain $\Omega_{HF}$, necessary for the offline stage of Proper Orthogonal Decomposition (POD) computed in $\Omega_{lf} \subset \Omega_{HF}$.

<figure><img src="image_6.png" /><figcaption>**FIGURE 5** | Snapshots of the water-air interface for the test case $H_\omega = 0.6$, over a bi-dimensional section of $\Omega_{HF}$. (a) Snapshot at $t = 0$ s. (b) Snapshot at $t = 1$ s. (c) Snapshot at $t = 2$ s. (d) Snapshot at $t = 3$ s. [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>
<table><caption>TABLE 1 | Mesh characteristics and relative errors with respect to the finest mesh $M_1$.</caption><thead><tr><th>Mesh</th><th>$h(m)$</th><th>$n_c$</th><th>$-\bar{\varepsilon}_{\text{mesh}}^u$</th><th>$-\bar{\varepsilon}_{\text{mesh}}^v$</th><th>$-\bar{\varepsilon}_{\text{mesh}}^p$</th><th>$-\bar{\varepsilon}_{\text{mesh}}^\alpha$</th><th>$-\bar{\varepsilon}_{\text{mesh}}^\beta$</th></tr></thead><tbody><tr><td>$M_4$</td><td>0.2</td><td>$1.8 \times 10^4$</td><td>1.7e-3</td><td>6.6e-4</td><td>3.9e-2</td><td>1.8e-4</td><td>2.1e-5</td></tr><tr><td>$M_3$</td><td>0.1</td><td>$8.5 \times 10^4$</td><td>1.0e-3</td><td>3.3e-4</td><td>1.5e-2</td><td>1.1e-4</td><td>1.0e-5</td></tr><tr><td>$M_2$</td><td>0.05</td><td>$4.4 \times 10^5$</td><td>5.2e-4</td><td>1.7e-4</td><td>8.6e-3</td><td>5.0e-5</td><td>3.5e-6</td></tr><tr><td>$M_1$</td><td>0.025</td><td>$2.5 \times 10^6$</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr></tbody></table>

Four different meshes {$M_i$}$_{i=1}^4$ are employed to discretize the domain $\Omega_{HF}$ spatially. In all cases, a hierarchical Cartesian (octree) mesh is utilized. Far from the air-water interface, the mesh is uniform in all three directions, that is $\Delta x = \Delta y = \Delta z = h$. Each mesh includes two levels of refinement near the air-water interface, with the initial position defined by $H_{\text{w}} = 0.6$ m and $H_{\text{s}} = 0$ m in Equation (41). The simulation time-step is set to $\Delta t = 0.002$ s. The characteristics of the different meshes are summarized in Table 1, where $n_c^1 > n_c^2 > n_c^3 > n_c^4$, represents the initial number of cells for mesh $M_i$. Since Adaptive Mesh Refinement (AMR) is applied, the number of mesh cells may slightly evolve over time as the air-water interface moves. The average relative errors

$$ \bar{\varepsilon}_{\text{mesh}}^{\varepsilon} = \int_{T} \frac{\sqrt{\int_{\Omega_{\text{HF}}} \|\xi_{M_i} - \xi_{M_1}\|_2^2} \mathrm{d}\mathbf{x}}{\sqrt{\int_{\Omega_{\text{HF}}} \|\xi_{M_1}\|_2^2} \mathrm{d}\mathbf{x}} \mathrm{d}t $$

for $i = \{2, 3, 4\}$ and $\xi = u, v, w, p$, and $\alpha$, where integrals over $\Omega_{HF}$ are discretized on a uniform Cartesian grid that fits the coarsest level of the coarsest grid ($M_4$), are also reported in this table.

For each quantity ($u, v, w, \alpha$, and $p$), a first-order behavior is observed.

In the subsequent analysis, we will use the mesh $M_2$ with $h = 0.05$ m that presents a good trade-off between computational cost and accuracy.

## 3.2 | Sensitivity of the Multi-Fidelity Model Resolution

Once the solution $U_h$ is obtained on the $M_2$ mesh, the POD basis functions {$\Phi_i$}$_{i=1}^{N_r}$ can be computed in an offline process, and the coupling model (25–32) can then be integrated in time. The final goal is to compute, at each time step, the POD coefficients {$a_i^*$}$_{i=1}^{N_r}$.

The Navier-Stokes Equations (25–27) are discretized in time using an implicit scheme. It is thus natural to impose the boundary conditions (28–31) at time $t^{(n+1)} = t^{(n)} + \Delta t$, and the minimization problem becomes a fully coupled non-linear problem, where all Equations (25–32) have to be solved at the same time, that is the solution $U = (u, v, w, \alpha)$ and POD coefficients $a^*$ must be solved simultaneously. This non-linear minimization can be expensive in the numerical viewpoint, and the goal is, thus, to compare the solution of the implicit coupling with a more simple explicit one, where boundary conditions (28–31) are provided in an explicit way, at time $t^{(n)}$. The minimization problem (32) (or (34)) for an explicit coupling reduces to a least square problem, where the $N_r$ POD coefficients {$a_i^*$}$_{i=1}^{N_r}$ are computed

from the known solution $U$ obtained at the previous time step. As already mentioned, the number of discrete points $N_p$ over the overlapping domain $\Omega_o$ should satisfy $N_p > N_r - M$. Then, the least square problem requires to compute $N_r$ unknowns using $N_p + M$ equations. In the algebraic point of view, the over-determined minimization problem (34) is approximated by $a^* = \text{argmin}_a (Aa - b)^2$, where $A \in \mathbb{R}^{(N_p+M) \times N_r}$ and $b \in \mathbb{R}^{(N_p+M)}$.

The first $N_p$ lines of $A$ and $b$ are built from the term

$$ \frac{1}{|\Omega_o|} \int_{\Omega_o} (\alpha_{\text{hf}} - \tilde{\alpha})^2 \mathrm{d}\mathbf{x} $$

while the last $M$ lines are built from the term

$$ \frac{\beta}{M} \sum_{i=1}^{M} (\hat{\alpha}_i - \tilde{\alpha}_i)^2 $$

Numerically, we have

$$ A = \begin{pmatrix} \phi_1^a(x_1) & \phi_2^a(x_1) & \cdots & \phi_{N_r}^a(x_1) \\ \vdots & \vdots & \ddots & \vdots \\ \phi_1^a(x_{N_p}) & \phi_2^a(x_{N_p}) & \cdots & \phi_{N_r}^a(x_{N_p}) \\ \phi_1^a(x_1^S) & \phi_2^a(x_1^S) & \cdots & \phi_{N_r}^a(x_1^S) \\ \vdots & \vdots & \ddots & \vdots \\ \phi_1^a(x_M^S) & \phi_2^a(x_M^S) & \cdots & \phi_{N_r}^a(x_M^S) \end{pmatrix} \quad \text{and} \quad b = \begin{pmatrix} \alpha_{\text{hf}}(x_1, t_j) - \bar{\alpha}(x_1) \\ \vdots \\ \alpha_{\text{hf}}(x_{N_p}, t_j) - \bar{\alpha}(x_{N_p}) \\ \hat{\alpha}(x_1^S, t_j) - \bar{\alpha}(x_1^S) \\ \vdots \\ \hat{\alpha}(x_M^S, t_j) - \bar{\alpha}(x_M^S) \end{pmatrix} \quad (42) $$

where the mean is subtracted from the high-fidelity solution to define $b$, which is a vector calculated for each $j$-th time instant. We recall that the interpolator $I_s$ must be used to compute the vector $b$ when the sensor locations do not coincide with the points on the low-fidelity mesh.

The POD coefficients are obtained by solving the associated normal equation $A^T A a^* = A^T b$.

Initially, to test the coupling algorithm, a trial without the POD model is conducted. Therefore, the boundary conditions (28–31) are provided by the reference high-fidelity solution $\tilde{U} = U_{\text{HF}}$, which is pre-computed on the large domain for the same physical configuration (41). This approach avoids truncation errors that arise from the POD expansion ($N_r \le N_s$), ensuring that only errors from the interpolation step are introduced into the model. Essentially, this method is equivalent to considering all POD modes constructed from snapshots computed at all times $t^{(k)}$, without any truncation, and sets a lower bound for the model error.

In the implicit version, conditions (28–31) are imposed using $\tilde{U}^{(n+1)} = U_{\text{HF}}^{(n+1)}$, while in the explicit version, the
conditions (28–31) are imposed from $\tilde{U}^{(n)} = U_{\text{HF}}^{(n)}$. In both versions, due to possible different meshes used in the large computational domain $\Omega_{\text{HF}}$ and in the small computational domain $\Omega_{\text{hf}}$, interpolations are performed, and we have $\tilde{U}^{(k)} = I_{\downarrow}(U_{\text{HF}}^{(k)})$, where $k = n$ for the explicit version, and $k = n + 1$ for the implicit version. The boundary conditions imposed for the small computation domain (28–31) are thus $U^{(k)} = I_{\uparrow} \circ I_{\downarrow}(U_{\text{HF}}^{(k)})$. Unfortunately, on the numerical viewpoint, we do not have $I_{\downarrow} = I_{\uparrow}^{-1}$, and numerical errors are, thus, necessarily introduced, leading to a non-zero error even if appropriate implicit boundary conditions are applied. Moreover, due to interpolations on the velocity fields, the mass conservation may be locally violated, possibly leading to some pressure oscillations.

For the following numerical example, a high-fidelity solution $U_{\text{HF}}$ (velocity, VOF, and pressure) is first computed on the entire large domain $\Omega_{\text{HF}}$ presented in Figure 4 up to $T = 3$ s, and snapshots are extracted at each computational time step $\Delta t$, chosen to be fixed in this study. A high-fidelity solution $U_{\text{hf}}$ is then computed in the small high-fidelity domain $\Omega_{\text{hf}}$, where boundary conditions $U = I_{\uparrow} \circ I_{\downarrow}(U_{\text{HF}})$ are applied on $\Gamma_{\text{hf}}$.

From now, and until the sensitivity analysis of the size of the domain $\Omega_{\text{hf}}$, we consider the case $R_6$ that will be presented in Section 3.4. No overlapping region is needed here.

Two cases are considered for this multi-fidelity coupling: a time step $\Delta t = 0.002$ s and a time step $\Delta t = 0.004$ s. For both time steps, the errors (39) (computed for the velocity, the volume fraction, and for the pressure) are comparable in the implicit and explicit cases; see Table 2. Due to the interpolations, the errors using implicit boundary conditions are not zero. Since implicit and explicit coupling give comparable errors, for sake of simplicity and efficiency, the explicit version with $\Delta t = 0.004$ s will be used in what follows. It should be noted that the error in the $z$ direction is larger compared to the other directions. This is primarily because the velocity component does not evolve in the $z$ direction due to the problem being predominantly two-dimensional, and the domain having limited spatial extension in the third direction. Consequently, the absolute errors in the $z$ direction are very small. However, it should be kept in mind that, since the presented errors are all relative errors, the ratio results in a higher overall value for $\bar{\epsilon}_{\text{hf}}^{w}$, and for the rest of the paper. Moreover, as expected, there are pressure oscillations, due to the interpolation process. However, despite these oscillations being noticeable, they are limited in magnitude, leaving the overall relative error very small, as demonstrated in Table 2. A possible strategy to reduce such pressure oscillations may be the implementation of a divergence-free interpolation step in the coupling algorithm.

For in-sample test, that is when same parameters in (41) are used in both $\Omega_{\text{HF}}$ and $\Omega_{\text{hf}}$, the error (37) computed with exact explicit boundary condition gives a lower bound for the coupling methodology. The exact boundary conditions can be equivalently computed by POD using $N_r = N_s = N_i$, with $N_r$ the number of POD modes, $N_s$ the number of snapshots, and $N_i$ the number of temporal iterations required by the high-fidelity solution. Extra errors will, thus, inevitably be added when considering reduced basis with $N_r \ll N_s \le N_i$.

### 3.3 | POD Sensitivity Analysis

The goal now is to perform sensitivity analysis of the number of retained POD modes $N_r$, and of the number of snapshots $N_s$. A numerical simulation is run in domain $\Omega_{\text{HF}}$ with $H_w = 0.60$ m and $H_s = 0$ m, up to $T = 3$ s. For the sensitivity of the number of POD modes $N_r$, snapshots of $U_{\text{HF}} = (u, v, w, \alpha)^T$ are collected at each time-step $t_j = j \Delta t$ with $\Delta t = 0.004$ s, and we have $N_s = N_i = 750$. The POD is computed from modified snapshots $U = U_{\text{HF}} - \bar{U}$, where $\bar{U}$ denotes the temporal average over the $N_s$ snapshots.

The POD eigenvalue spectrum $\{\lambda_i/\lambda_1\}_{i=1}^{N_s}$ and relative information content $\text{RIC}(N_r) = \sum_{i=1}^{N_r} \lambda_i / \sum_{i=1}^{N_s} \lambda_i$ are plotted in Figure 6.

The relative projection errors for $\xi = u, v, w$ and $\alpha$ over the domain $\Omega_{\text{HF}}$

$$ \bar{\epsilon}_{\text{proj}}^{\xi}(N_r) = \int_T \frac{\sum_{j=1}^{N_r} \| I_{\downarrow}(\xi_{\text{HF}}(t_j)) - (\bar{\xi} + \sum_{k=1}^{N_r} a_k^{\text{proj}}(t_j) \Phi_k^{\xi}) \|_{\Omega_{\text{HF}}}}{\sum_{j=1}^{N_r} \| I_{\downarrow}(\xi_{\text{HF}}(t_j)) \|_{\Omega_{\text{HF}}}} \, dt \quad (43) $$

are plotted in Figure 7, where the temporal coefficients $a_k^{\text{proj}}(t_j)$ correspond to the POD eigenvectors. These errors traduce the filtering effect of the POD truncation, and quantify, thus, the error due to the truncated modes.

The trend of the errors plotted in Figure 7 follows the trend of the POD eigenvalues in Figure 6: the errors are decreasing with the number of POD modes. The trend is less pronounced for the $w$ velocity component due to the previously mentioned bidimensional nature of the problem.

In what follows, a sensitivity analysis of the number of retained POD modes, $N_r$, and the number of snapshots, $N_s$, is performed on the hybrid multi-fidelity model (25–27). The errors (37) and (38) are computed for velocity and VOF.

<table><caption>TABLE 2 | Relative errors of the multi-fidelity approach with two time steps, for explicit and implicit coupling.</caption><thead><tr><th>Domain</th><th>Case</th><th>ε<sup>-u</sup><sub>hf</sub></th><th>ε<sup>-v</sup><sub>hf</sub></th><th>ε<sup>-w</sup><sub>hf</sub></th><th>ε<sup>-α</sup><sub>hf</sub></th><th>ε<sup>-p</sup><sub>hf</sub></th></tr></thead><tbody><tr><td rowspan="2">Ω<sub>hf</sub></td><td>Implicit</td><td>dt = 0.002 s</td><td>2.60e-5</td><td>1.69e-5</td><td>1.89e-3</td><td>5.78e-6</td><td>2.06e-6</td></tr><tr><td>Explicit</td><td>dt = 0.004 s</td><td>7.15e-5</td><td>5.13e-5</td><td>3.19e-3</td><td>1.49e-5</td><td>4.89e-6</td></tr><tr><td></td><td>Implicit</td><td>dt = 0.002 s</td><td>2.98e-5</td><td>1.91e-5</td><td>1.94e-3</td><td>7.58e-6</td><td>2.37e-6</td></tr><tr><td></td><td>Explicit</td><td>dt = 0.004 s</td><td>8.67e-5</td><td>5.93e-5</td><td>3.28e-3</td><td>2.23e-5</td><td>6.00e-6</td></tr></tbody></table>
<figure><img src="image_7.png" /><figcaption>FIGURE 6 | Eigenvalue spectrum (left) and RIC (right) with respect to the number of POD modes $N_r$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIGURE 7 | Influence of the number of POD modes $N_r$ on the projection error $\epsilon_{proj}$.</figcaption></figure>

### 3.3.1 | Sensitivity of the Number of Modes $N_r$

The sensitivity of the number of modes $N_r$ is studied keeping constant the number of snapshots $N_s = 750$ over the entire domain $\Omega_{HF}$. We recall that the goal of the hybrid fidelity coupling is twofold. The POD reduced order solution $\tilde{U} = \bar{U} + \Phi a^*$ with the optimized coefficients $a^*$ is used to give appropriate boundary conditions for the small high-fidelity domain $\Omega_{hf}$, and also to propagate information through the large domain, $\Omega_{lf} \setminus \Omega_o = \Omega_{HF} \setminus \Omega_{hf}$.

From now, and until the sensitivity analysis of the size of the overlapping domain $\Omega_o$ that will be presented in Section 3.5, we consider $\Omega_o = \Omega_{hf}$.

Both errors (39) and (40) are reported in Table 3. One of the most important variables for the future wave energy converter

application is the VOF $\alpha$, and Table 3 shows that starting from $N_r = 20$, the VOF error is almost constant in the domain $\Omega_{hf}$. Although theory suggests that increasing the number of modes would reduce the error, there is likely a filtering effect caused by the interpolator $I_1$ used to impose boundary conditions. Since the interpolation step is essential for the coupling methodology as presented, and given the limitations of the commercial software's interpolation options, this filtering effect is not ultimately detrimental. It allows us to use a small $N_r$ while still achieving good results. Since all the other errors are still decreasing with $N_r$, we arbitrarily chose a slightly larger value, $N_r = 30$.

A comparison of the temporal evolutions of the first five POD coefficients $a_k^{\text{proj}}(t)$ obtained by projection (eigenvectors of the POD) and by the coupling strategy $a_k^*(t)$, is given in Figure 8. Good agreements are observed, validating, thus, the reduced order model obtained with the Galerkin-free coupling strategy.

<table><caption>TABLE 3 | Relative errors over the high-fidelity domain $\Omega_{hf}$ and over the remaining low-fidelity domain $\Omega_{lf} \setminus \Omega_o$.</caption><thead><tr><th>Domain</th><th>Case</th><th>&epsilon;&#773;<sub>u</sub></th><th>&epsilon;&#773;<sub>v</sub></th><th>&epsilon;&#773;<sub>w</sub></th><th>&epsilon;&#773;<sub>&alpha;</sub></th></tr></thead><tbody><tr><td rowspan="6">&Omega;<sub>lf</sub> \ &Omega;<sub>o</sub></td><td>N<sub>r</sub> = 5</td><td>7.28e-4</td><td>6.53e-4</td><td>1.25e-3</td><td>1.10e-4</td></tr><tr><td>N<sub>r</sub> = 10</td><td>2.57e-4</td><td>1.88e-4</td><td>9.98e-4</td><td>4.75e-5</td></tr><tr><td>N<sub>r</sub> = 20</td><td>1.17e-4</td><td>8.09e-5</td><td>7.98e-4</td><td>2.60e-5</td></tr><tr><td>N<sub>r</sub> = 30</td><td>7.63e-5</td><td>5.35e-5</td><td>7.09e-4</td><td>1.86e-5</td></tr><tr><td>N<sub>r</sub> = 50</td><td>6.01e-5</td><td>4.61e-5</td><td>6.33e-4</td><td>1.28e-5</td></tr><tr><td>N<sub>r</sub> = 100</td><td>5.71e-5</td><td>4.68e-5</td><td>5.95e-4</td><td>1.11e-5</td></tr><tr><td rowspan="6">&Omega;<sub>hf</sub></td><td>N<sub>r</sub> = 5</td><td>1.36e-3</td><td>1.10e-3</td><td>1.06e-1</td><td>2.61e-4</td></tr><tr><td>N<sub>r</sub> = 10</td><td>4.42e-4</td><td>3.21e-4</td><td>2.11e-2</td><td>6.55e-5</td></tr><tr><td>N<sub>r</sub> = 20</td><td>3.34e-4</td><td>2.27e-4</td><td>2.91e-2</td><td>2.58e-5</td></tr><tr><td>N<sub>r</sub> = 30</td><td>1.91e-4</td><td>1.32e-4</td><td>1.08e-2</td><td>2.63e-5</td></tr><tr><td>N<sub>r</sub> = 50</td><td>1.49e-4</td><td>1.08e-4</td><td>7.81e-3</td><td>2.65e-5</td></tr><tr><td>N<sub>r</sub> = 100</td><td>1.16e-4</td><td>8.40e-5</td><td>4.86e-3</td><td>2.73e-5</td></tr></tbody></table>
### 3.3.2 | Sensitivity of the Number of Snapshots $N_s$

The number of POD basis functions is $N_r = 30$. The projection errors with respect to $N_s$ are shown in Table 4. For each variable, a first order is observed: the error is divided by $n$, using $n$ times more snapshots. We arbitrarily choose $N_s = 250$ in what follows.

## 3.4 | Sensitivity of the Size of High-Fidelity Domain $\omega_{hf}$

An important choice concerns the size of the small high-fidelity domain $\Omega_{hf}$. We still consider here $\Omega_o = \Omega_{hf}$. Each configuration under consideration is noted $R_D$, and is parameterized with the ratio $D = |\Omega_{hf}|/|\Omega_{hf}|$. The definition of $R_D$ for the sensitivity analysis is constrained by the simulation setup. The lower bound, indicating

the smallest size of $\Omega_{hf}$ (or higher $D$), is determined by the size of the body being simulated. In the context of placing a floating circular body, $\Omega_{hf}$ needs to be sufficiently large to fully encompass the body in the high-fidelity domain, including the overset mesh cells around it. The upper bound is chosen in order to have a minimum of reduction in the computational cost. We consider here two extreme sizes, namely a large one with $R_3$, and a small one with $R_8$. Another size is also considered, $R_D = 6$, which is a trade-off that significantly reduces computational cost, while ensuring the body remains completely within $\Omega_{hf}$, and is used in all the simulations previously shown. The three cases are outlined in Figure 9. Table 5 shows the ratio of the number of mesh cells at initialization, $R_{n_{cells}} = \frac{\Omega_{hf}}{\Omega_{hf}}$. The relation is not linear, but depends on the form and position of $\Omega_{hf}$, because the majority of the cells are concentrated around the free surface, where there is a mesh refinement.

The errors (39) and (40) are displayed in Table 6. We can observe that the errors obtained for the $R_3$ and $R_6$ cases are almost the same, while the error for the $R_8$ case is larger. A good trade-off between computational time and accuracy seems to be the $R_6$ case. This case has been chosen in the previous sections and will continue to be used in what follows.

<figure><img src="image_9.png" /><figcaption>FIGURE 8 | Comparison of the POD temporal coefficients $a^{proj}(t)$ obtained by projection (eigenvectors of the POD) and by the coupling strategy, $a^*(t)$. [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

<table><caption>TABLE 4 | Sensitivity analysis on the number of snapshots $N_s$.</caption><thead><tr><th rowspan="2">N<sub>s</sub></th><th colspan="2">Multiple</th><th rowspan="2">RIC (%)</th><th rowspan="2">ε<sup>-u</sup><sub>proj</sub></th><th rowspan="2">ε<sup>-v</sup><sub>proj</sub></th><th rowspan="2">ε<sup>-w</sup><sub>proj</sub></th><th rowspan="2">ε<sup>-α</sup><sub>proj</sub></th></tr><tr><th>of Δt</th><th>Δt</th></tr></thead><tbody><tr><td>50</td><td>15</td><td>15</td><td>99.80</td><td>6.4e-4</td><td>2.6e-4</td><td>7.3e-3</td><td>2.7e-4</td></tr><tr><td>75</td><td>10</td><td>10</td><td>99.78</td><td>4.6e-4</td><td>2.0e-4</td><td>5.6e-3</td><td>1.9e-4</td></tr><tr><td>150</td><td>5</td><td>5</td><td>99.77</td><td>2.3e-4</td><td>1.1e-4</td><td>3.1e-3</td><td>9.7e-5</td></tr><tr><td>250</td><td>3</td><td>3</td><td>99.77</td><td>1.4e-4</td><td>6.4e-5</td><td>1.9e-3</td><td>5.8e-5</td></tr><tr><td>750</td><td>1</td><td>1</td><td>99.77</td><td>4.7e-5</td><td>2.1e-5</td><td>6.4e-4</td><td>1.9e-5</td></tr></tbody></table>

## 3.5 | Sensitivity of the Size of the Overlapping Domain $\omega_o$

In the previous analysis, the overlapping region was chosen to be $\Omega_o = \Omega_{hf}$, which corresponds hereafter to the reference case 0. Three other frame-shaped domains $\Omega_o$ are also considered and are presented in Figure 10.

If no sensors are considered ($M = 0$), the minimization problem for $a^*$ is well-posed if the number of mesh cells $N_p$ included in

<table><caption>TABLE 5 | Number of mesh cells ratio $R_{n_{cells}}$ for the three test cases $R_3$, $R_6$, and $R_8$.</caption><thead><tr><th>Case</th><th>R<sub>n<sub>cells</sub></sub></th></tr></thead><tbody><tr><td>R<sub>3</sub></td><td>1.69</td></tr><tr><td>R<sub>6</sub></td><td>3.38</td></tr><tr><td>R<sub>8</sub></td><td>4.50</td></tr></tbody></table>

<figure><img src="image_10.png" /><figcaption>FIGURE 9 | Test cases for varying size of $\Omega_{hf}$. [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>
$\Omega_o$ is $N_p \ge N_r$. Each configuration considered here satisfies this requirement.

The average relative errors computed in the small high-fidelity domain $\Omega_{hf}$ and in the low-fidelity domain $\Omega_{lf} \setminus \Omega_o$ are reported in Table 7. The errors obtained in case (c) are quite close to those obtained in case 0 (used until now). The errors obtained in cases (a) and (b) are one order of magnitude higher. In this problem, the primary information is concentrated along the x-direction, with the most critical boundary conditions for $\Omega_{hf}$ being the lateral ones (left and right), where appropriate conditions, particularly for the water level, must be imposed numerically using the VOF $\alpha$. Consequently, an overlapping region shaped like that in case (c), with a greater extent in the longitudinal direction than in the vertical, provides better results, and will be used in what follows. Changing the size and location of the overlapping region could result in different values of the temporal coefficients, and thus, different overall outcomes.

<table><caption>TABLE 6 | Errors (39) and (40) for test cases $R_3$, $R_6$, and $R_8$.</caption><thead><tr><th>Domain</th><th>Case</th><th>-ε<sup>u</sup></th><th>-ε<sup>v</sup></th><th>-ε<sup>w</sup></th><th>-ε<sup>-α</sup></th></tr></thead><tbody><tr><td rowspan="3">Ω<sub>lf</sub> \ Ω<sub>o</sub></td><td>R<sub>3</sub></td><td>6.77e-5</td><td>3.92e-5</td><td>7.38e-4</td><td>1.42e-5</td></tr><tr><td>R<sub>6</sub></td><td>7.63e-5</td><td>5.35e-5</td><td>7.09e-4</td><td>1.86e-5</td></tr><tr><td>R<sub>8</sub></td><td>1.34e-4</td><td>1.20e-4</td><td>7.86e-4</td><td>2.52e-5</td></tr><tr><td rowspan="3">Ω<sub>hf</sub></td><td>R<sub>3</sub></td><td>1.44e-4</td><td>1.41e-4</td><td>9.69e-3</td><td>2.08e-5</td></tr><tr><td>R<sub>6</sub></td><td>1.91e-4</td><td>1.32e-4</td><td>1.08e-2</td><td>2.63e-5</td></tr><tr><td>R<sub>8</sub></td><td>6.31e-4</td><td>3.09e-4</td><td>2.27e-2</td><td>4.88e-5</td></tr></tbody></table>

## 3.6 | Sensitivity of the Sampling for Out-Of-Sample Prediction

The final configurations resulting from the previous sensitivity analysis are as follows:

* The $M_2$ mesh is used to discretize the large high-fidelity computational domain $\Omega_{HF}$.

* The time step is $\Delta t = 0.004$ s.

* An explicit coupling is used to exchange information between high- and low-fidelity models.

* The number of snapshots is $N_s = 250$ (extracted uniformly with $\Delta t_s = 0.012$ s over $T = 3$ s).

* The number of retained POD modes is $N_r = 30$.

* The small high-fidelity domain $\Omega_{hf}$ corresponds to case $R_6$.

* The overlapping domain $\Omega_o$ corresponds to case (c).

<table><caption>TABLE 7 | Relative errors with respect to the overlapping domain.</caption><thead><tr><th>Domain</th><th>Case</th><th>-ε<sup>u</sup></th><th>-ε<sup>v</sup></th><th>-ε<sup>w</sup></th><th>-ε<sup>-α</sup></th></tr></thead><tbody><tr><td rowspan="4">Ω<sub>lf</sub> \ Ω<sub>o</sub></td><td>0</td><td>7.61e-5</td><td>5.33e-5</td><td>7.22e-4</td><td>1.86e-5</td></tr><tr><td>a</td><td>2.23e-3</td><td>1.98e-3</td><td>4.45e-3</td><td>2.92e-4</td></tr><tr><td>b</td><td>1.81e-3</td><td>1.80e-3</td><td>2.47e-3</td><td>2.59e-4</td></tr><tr><td>c</td><td>2.30e-4</td><td>2.06e-4</td><td>8.05e-4</td><td>2.83e-5</td></tr><tr><td rowspan="4">Ω<sub>hf</sub></td><td>0</td><td>1.85e-4</td><td>1.29e-4</td><td>1.14e-2</td><td>2.61e-5</td></tr><tr><td>a</td><td>3.50e-3</td><td>3.63e-3</td><td>4.28e-1</td><td>3.67e-4</td></tr><tr><td>b</td><td>2.48e-3</td><td>2.32e-3</td><td>1.39e-1</td><td>5.12e-4</td></tr><tr><td>c</td><td>3.09e-4</td><td>2.49e-4</td><td>1.17e-2</td><td>4.73e-5</td></tr></tbody></table>

<figure><img src="image_11.png" /><figcaption>FIGURE 10 | Overlapping domains $\Omega_o$ tested.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**FIGURE 11** | Quantitative (left) and graphical (right) definition of the variable parameter, $H_w$, for the out-of-sample tests. Case 0 is the in-sample test, taken as reference. Case E groups the snapshots from all the other cases, except Case 0. [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

In the previous sections, the data collected from simulations with initial conditions using $H_w = 0.60$ m are used to build the POD basis functions for the coupled model that reproduces the simulation for the same initial conditions. In other words, we are performing an in-sample test, also called a reproduction problem. In this section, data from simulations based on different initial conditions $H_w$ are collected and used in the training phase, and the coupling algorithm is evaluated for the simulation initialized with $H_w = 0.60$ m (case 0), which is not included in the database. Several training sets are considered (see Figure 11), combining different $H_w$. The chosen $H_w$ are progressively farther from the target initial conditions (cases A, B, C, and D), or are all collected in a large dataset (case E). Note that we use $N_s = 250$ snapshots for each simulation, and, for instance, the database for case E is thus composed by $8 \times 250$ snapshots. However, for each case (0, A, B, C, D, and E), $N_r = 30$ POD modes are considered.

Table 8 shows that the solution for the baseline configuration (case 0) can be well approximated using training sets that do not contain the in-sample snapshots. In particular, datasets from initial conditions closer to the target (case A) give good results. The relative errors increase when the training points are farther from the baseline configuration. The extensive dataset (case E) yields satisfactory results, albeit with higher errors compared to case A. This outcome arises due to an overflow of information: case E incorporates not only critical data from case A, but also less pertinent information from cases B-D, which can diminish the overall performance. However, the benefit of employing a dataset like case E lies in its enhanced robustness across a broad spectrum of $H_w$ values, compared to cases A-D. This robustness stems from the system being trained on numerous $H_w$ values, thereby making case E highly versatile. Better sampling strategies can be considered, but it is out of the scope of this paper.

## 4 | Wave Energy Converter

The goal is to numerically model a class of wave energy converters, where the energy harvesting devices are floating bodies with rigid motions driven by sea waves. Floating wave energy converters are typically deployed in deep waters, ideally located far from the shoreline. These devices are equipped with a mooring system

<table><caption>**TABLE 8** | Relative errors on the high-fidelity domain $\Omega_{hf}$ and on the remaining low-fidelity domain $\Omega_{lf} \setminus \Omega_o$ for the different sampling test cases.</caption><thead><tr><th>Domain</th><th>Case</th><th>ε̄<sub>u</sub></th><th>ε̄<sub>v</sub></th><th>ε̄<sub>w</sub></th><th>ε̄<sub>α</sub></th></tr></thead><tbody><tr><td rowspan="6">Ω<sub>lf</sub> \ Ω<sub>o</sub></td><td>0</td><td>6.29e-5</td><td>4.23e-5</td><td>7.09e-4</td><td>1.74e-5</td></tr><tr><td>A</td><td>1.09e-4</td><td>8.57e-5</td><td>1.27e-3</td><td>2.41e-5</td></tr><tr><td>B</td><td>2.68e-4</td><td>2.32e-4</td><td>1.43e-3</td><td>4.84e-5</td></tr><tr><td>C</td><td>3.54e-4</td><td>2.82e-4</td><td>1.47e-3</td><td>6.75e-5</td></tr><tr><td>D</td><td>6.76e-4</td><td>5.56e-4</td><td>1.57e-3</td><td>1.18e-4</td></tr><tr><td>E</td><td>3.21e-4</td><td>2.68e-4</td><td>1.30e-3</td><td>5.83e-5</td></tr><tr><td rowspan="6">Ω<sub>hf</sub></td><td>0</td><td>1.72e-4</td><td>1.19e-4</td><td>1.26e-2</td><td>1.37e-5</td></tr><tr><td>A</td><td>3.43e-4</td><td>2.14e-4</td><td>2.29e-2</td><td>1.95e-5</td></tr><tr><td>B</td><td>3.64e-4</td><td>2.48e-4</td><td>2.11e-2</td><td>5.14e-5</td></tr><tr><td>C</td><td>4.50e-4</td><td>2.98e-4</td><td>1.66e-2</td><td>5.97e-5</td></tr><tr><td>D</td><td>6.18e-4</td><td>4.94e-4</td><td>2.26e-2</td><td>2.02e-4</td></tr><tr><td>E</td><td>4.58e-4</td><td>2.96e-4</td><td>2.04e-2</td><td>7.12e-5</td></tr></tbody></table>

and harness the oscillatory motion generated by the waves to produce energy [58]. They can exhibit varying degrees of freedom, with notable prototypes including point absorbers, overtopping devices, surface attenuators, and rotating mass devices.

We consider here a sphere floating on sea waves, the latter being imposed by $M > 0$ sensors in the coupling strategy. The floating body is defined as a point absorber, that is the only degree of freedom is in the y direction (aligned with gravity). The objective is to predict the wave behavior and the motion of the floating body subjected to an imposed wave. The wave is defined with some characteristics $W0$ (period, height, and wavelength) defined in Table 9. The sketch of the domain decomposition is depicted in Figure 12. The boundary $\Gamma_{hf}$ now has an extra part, orthogonal to the z direction (see top figure in Figure 12), and the green points represent the $M$ sensors where information on the incoming wave is measured. Such sensor area is defined in this study to be symmetrically placed with respect to the incoming wave in the vertical direction, and as compact as possible, while still encompassing a few cells, longitudinally. In this test, $\beta$ is set to unity,
equalizing the influence of both the floating body and the wave. Alternatively, adjusting the value of $\beta$ allows for varying the influence of the waves on the converter, or vice versa.

The initial computational domain, noted $\Omega_{sim}$, expands to 85 × 45 × 10 m. The large size of this domain is necessary to take into account the water depth, the multiple wavelengths, and artificial forcing zones that avoid spurious reflections. A symmetry condition is used to reduce the computational burden of the FOM simulation, and is also kept in the coupled model; see Figure 12.

However, without loss of generality, we consider a smaller domain $\Omega_{HF} \subset \Omega_{sim}$ to reduce the memory footprint and CPU costs for POD basis functions. Even though it is necessary to compute the actual solution on $\Omega_{sim}$, we only consider its trace on $\Omega_{HF}$, which is a subset of the larger domain, and where no further computations are required. Similar to the previous test case, the domain $\Omega_{HF}$ is 10 × 6 × 0.8 m. This allows us to consider more than one wavelength in the domain. However, the size of $\Omega_{HF}$ in the z direction is twice the one used in the sensitivity analysis section. This is necessary to numerically model three-dimensional effects near the floating body. Moreover, the domain has only one symmetry plane, passing through the middle of the sphere. The dimension of $\Omega_{hf}$ also slightly changes. Even though the x and y dimensions still correspond to the $R_6$ case, there is a modification in the z direction, with the creation of the new boundary (see top of Figure 12). Likewise, the overlapping region conserves the same size in the x and y directions, corresponding to test case (c) presented in the dedicated section, but a new area is added in the z direction. With the new dimensions, the volume ratio between $\Omega_{HF}$ and $\Omega_{hf}$ becomes 8, yielding a bigger reduction with respect to the test case used in

the sensitivity analysis. The domain $\Omega_{hf}$ is 4 × 2.5 × 0.6 m, which represents 0.015% of the volume of the real simulation domain ($\frac{|\Omega_{sim}|}{|\Omega_{hf}|} = 6375$).

Numerically, $\Omega_{HF}$ is composed of 413457 cells, $\Omega_{hf}$ is composed of 53703 cells, and $\Omega_o$ is composed of 38493 cells. For simplicity reasons, we use $M = 3482$ sensors covering a domain of dimensions 0.2 × 2 × 0.8 m; see Figure 14.

An overset mesh is used for the sphere in the whole domain $\Omega_{sim}$, and in the small domain $\Omega_{hf}$. The POD basis functions are only computed on $\Omega_{hf} = \Omega_o \cup (\Omega_{HF} \setminus \Omega_{hf})$.

We denote by IS the In-Sample, or reconstruction case, where the POD basis is computed with a database from the baseline condition $W0$, and by OOS the Out-Of-Sample case where the POD basis is computed with a database from conditions $W1$ and $W2$, and the solution is predicted for wave $W0$. Waves $W0$, $W1$ and $W2$ are defined in Table 9.

The number of POD modes is always $N_r = 30$. Each simulation lasts $T = 6$ s, allowing an incoming wave to impact the floater. The snapshots are collected at $\Delta t = 0.012$ s, for a total of $N_s = 500$ for each simulation. These settings are derived from the previous sensitivity analysis to ensure a significant reduction, while retaining enough information in the reduced model to accurately simulate the wave. Altering the parameters would likely change the outcome. For instance, reducing the number of modes would decrease the quality of the results, leading to a poor model with potential spurious, non-physical effects. In an out-of-sample test, this issue would be exacerbated if the training set's characteristics were very different from the tested wave. Therefore, despite the masking effect of the interpolation step, which limits quality improvement with more modes, there is still a minimum number of modes that must remain physically acceptable.

The computed errors are reported in Table 10. The errors obtained in the Out-Of-Sample (OOS) case are slightly larger than in the In-Sample (IS) case, but they are comparable. In this test case involving a spherical floating body, the third direction gains importance with respect to the previous test case, thus aligning

<table><caption>TABLE 9 | Wave characteristics.</caption><thead><tr><th>Wave</th><th>Wave period (s)</th><th>Wave height (m)</th><th>Wavelength (m)</th></tr></thead><tbody><tr><td>W0</td><td>2</td><td>0.2</td><td>6.25</td></tr><tr><td>W1</td><td>1.95</td><td>0.18</td><td>5.94</td></tr><tr><td>W2</td><td>2.05</td><td>0.22</td><td>6.56</td></tr></tbody></table>

<figure><img src="image_13.png" /><figcaption>FIGURE 12 | Domain decomposition for the sphere simulation. [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>
<table><caption>**TABLE 10** | Errors for In-Sample (IS) and Out-Of-Sample (OOS) tests.</caption><thead><tr><th>Domain</th><th>Case</th><th>ε̄<sub>u</sub></th><th>ε̄<sub>v</sub></th><th>ε̄<sub>w</sub></th><th>ε̄<sub>α</sub></th></tr></thead><tbody><tr><td rowspan="2">Ω<sub>lf</sub> \ Ω<sub>o</sub></td><td>IS</td><td>2.40e-4</td><td>1.67e-4</td><td>4.35e-4</td><td>3.02e-5</td></tr><tr><td>OOS</td><td>4.74e-4</td><td>3.39e-4</td><td>1.09e-3</td><td>5.23e-5</td></tr><tr><td rowspan="2">Ω<sub>hf</sub></td><td>IS</td><td>4.35e-4</td><td>3.08e-4</td><td>1.59e-3</td><td>4.38e-5</td></tr><tr><td>OOS</td><td>4.80e-4</td><td>3.87e-4</td><td>1.62e-3</td><td>7.77e-5</td></tr></tbody></table>

<figure><img src="image_14.png" /><figcaption>**FIGURE 13** | Evolution of the force acting on the body and the heave motion of the body. [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

the relative error for $w$ closer to that of the velocity variables in the other directions. However, due to the unidirectional nature of the waves, this influence remains minimal, and the error in $w$ continues to differ from the errors in other directions. The evolution of the vertical translation of the body and the vertical force acting on it are plotted in Figure 13. For the IS case, the temporal evolution of the force and the position of the body is very close to the reference one (HF), with a superposition of the curves, despite a loss in amplitude. In the OOS scenario, there is a slightly larger difference between the predicted solution and the reference. Indeed, the loss in amplitude is slightly higher than in the IS case, and, in addition, there is a phase shift of the curves; yet overall performance remains satisfactory. This small out-of-sample test demonstrates the potential for achieving accurate predictions with significant CPU savings. Future directions for this study include expanding the training set with additional waves and employing diverse clustering techniques, to create an extensive and robust training set while maintaining manageable sizes.

In terms of computational time, directly defining a speedup is not feasible since the original CFD simulation on $\Omega_{sim}$ and the coupled model are executed on different computers. However, certain observations can be made. Assessing the complete offline cost is challenging because the pre-processing involves various steps that depend on the specific test, computers, and codes used. However, excluding the creation of the snapshots,

<figure><img src="image_15.png" /><figcaption>**FIGURE 14** | Snapshot of the air-water interface (from the VOF variable $\alpha$) with the spherical point absorber, obtained with the multi-fidelity model in the OOS test case. The green domain represents $\Omega_{HF}$, the violet domains represent the total support of the minimization problem, given by both the overlapping domain (around the sphere) and the sensor locations, here taken over a region close to the inlet boundary (on the left). Only one half of the domain in the $z$ direction is actually computed, the visualization in the whole domain is obtained by symmetry. [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

processing each snapshot takes only a few seconds, and defining the basis functions takes a few minutes, on one processor. Therefore, depending on the simulation and the number of snapshots required, this part of the offline cost amounts to no more than a few hours. It is important to note that this is a one-time cost. Regarding the online cost, comparing the simulation time for the high-fidelity solvers offers a more direct indication of the advantages provided by the coupled methodology. A single high-fidelity simulation on $\Omega_{sim}$ with the wave characteristics outlined in Table 9 requires approximately 10 hours and 30 minutes using 48 processors. In contrast, the coupled model produces results in less than 4 hours using 6 processors. Although this speedup is significant, it does not reflect the storage gain ratio, as the computational bottleneck, implementing the DFBI model and the overset mesh for the moving body, is present in both high-fidelity simulations, $\Omega_{sim}$ and $\Omega_{hf}$. Furthermore, despite the necessary expense of the offline training phase, which involves running costly simulations for a few parameter values (here $W1$ and $W2$), it enables the generation of a wide array of new solutions for numerous other parameter values (with $W0$ representing just one example, and the most distant in terms of parameter values), at the reduced cost afforded by the coupled model.

Finally, a snapshot of the air-water interface (from the VOF variable $\alpha$) with the spherical point absorber, obtained with the multi-fidelity model in the OOS test case, is given in Figure 14. The air-water interface is continuous across the whole domain, no spurious effects being observed near the overlapping zone. The absence of reflections at the boundaries is most likely due to the stabilizing models implemented in the commercial software, allowing to impose the solution using only the boundaries,
without the need for buffer zones. Additionally, the precautions required for the Ωsim domain, which necessitates a very large simulation area for a correct definition of the wave, are not needed in Ωhf because only a multi-phase flow model is implemented, and not a wave model.

5 | Conclusion

The primary objective of this study was to efficiently model the behavior of Wave Energy Converters (WECs) when exposed to sea waves. This problem is intricate, as it involves the interactions between a floating structure and the fluid corresponding to actual sea waves. In this study, we considered a bi-fluid model for air and water, wherein the interface is embedded and modelled using the Volume Of Fluid (VOF) method. The Full Order Model (FOM), based on the Navier-Stokes equations with VOF and fluid-structure interactions, was solved using the commercial software StarCCM+. While the solution of the Full Order Model is accurate, it becomes prohibitively costly when applied to large-scale WEC problems. Indeed, modelling actual sea waves necessitates a significantly large computational domain to accurately generate and potentially dampen waves, particularly to prevent reflections.

In light of these challenges, Reduced Order Models (ROMs) are essential to save computational resources. In this study, we employed Proper Orthogonal Decomposition (POD) as the data compression reduction technique. However, classical POD projection-based (Petrov-)Galerkin ROM for bi-fluid flows poses efficiency issues with model order reduction. Its complexity, manifested in the form of a fourth-order tensor in the leading dynamical system, presents significant hurdles. To overcome this difficulty, hyper-reduction techniques can be used but are challenging or impossible to apply with commercial codes where the source code is not directly available.

To address these challenges, we have proposed a non-intrusive approach that does not require access to numerical residuals or mathematical models. We have developed a domain decomposition Galerkin-free approach coupled with sensor information to impose incoming waves. Our multi-fidelity model integrates POD ROM and FOM spatially, with FOM utilized near obstacles where viscous effects and non-linearities prevail, while a POD ansatz is employed for smoother flow regions.

Our coupling methodology simplifies to a least square minimization problem, making it easily applicable to simulations of moving bodies, including wave energy conversion.

We have began by considering a simple prototype, a point absorber sphere, but our methodology can easily be extended to more complex and full-scale WECs with multiple degrees of freedom. Furthermore, without altering the coupling strategy, additional details such as mooring or control systems can be incorporated for a more realistic simulation.

Parameter tuning, such as adjusting domain sizes or the number of POD modes, is problem-specific and can be tailored to strike a balance between accuracy and speed-ups. The accuracy is

improved considering a large number of POD modes with a large high-fidelity domain. Conversely, the speed-up can be increased considering a small high-fidelity domain coupled with few POD modes.

The approach has been validated for in-sample cases, where the POD subspace was trained using appropriate wave parameters. This represents a reproduction scenario, as we previously simulated this problem to obtain the POD modes. More notably, our method is also able to predict solutions for unseen parameters in out-of-sample simulations. The results maintain accuracy while significantly reducing CPU costs. This is achieved by leveraging a small high-fidelity domain, which represents only a fraction of the original domain used in training the POD modes. The efficiency of our hybrid fidelity model enables intensive simulations, making it suitable for optimization tasks and facilitating multiple-query simulations.

Our ongoing work includes addressing the optimal sensor location problem to further improve results. Additionally, we are exploring clustering techniques on snapshots to enhance efficiency. This involves reducing the POD offline stage while simultaneously improving the online stage by decreasing the number of POD modes for each cluster. This strategy can be viewed as a method for approximating a non-linear problem with piece-wise linear approaches.

Our ultimate goal is to model and optimize farms comprising multiple WECs. We envision achieving this through an algorithm that integrates multiple POD/FOM patches, augmenting the size of the minimization problem while retaining simplicity as a least-square problem.

Acknowledgments

The authors acknowledge on-premise HPC cluster resources made available by Optimad Srl for conducting the research reported in this paper. The authors acknowledge the support by the European Union's Horizon 2020 Research and Innovation Programme under the Marie Skłodowska-Curie Actions, grant agreement 872442 (ARIA).

Conflicts of Interest

The authors declare no conflicts of interest.

Data Availability Statement

The data that support the findings of this study are available from the corresponding author upon reasonable request.

References

1. G. S. Payne, J. R. Taylor, T. Bruce, and P. Parkin, "Assessment of Boundary-Element Method for Modelling a Free-Floating Sloped Wave Energy Device. Part 2: Experimental Validation," *Ocean Engineering* 35, no. 3 (2008): 342–357, https://doi.org/10.1016/j.oceaneng.2007.10.008.

2. A. Rahimi, S. Rezaei, J. Parvizian, et al., "Numerical and Experimental Study of the Hydrodynamic Coefficients and Power Absorption of a Two-Body Point Absorber Wave Energy Converter," *Renewable Energy* 201 (2022): 181–193, https://doi.org/10.1016/j.renene.2022.10.103.
3. G. Giorgi and J. V. Ringwood, "Nonlinear Froude-Krylov and Viscous Drag Representations for Wave Energy Converters in the Computation/Fidelity Continuum," *Ocean Engineering* 141 (2017): 164-175, https://doi.org/10.1016/j.oceaneng.2017.06.030.

4. A. Engsig-Karup, C. Eskilsson, and D. Bigoni, "A Stabilised Nodal Spectral Element Method for Fully Nonlinear Water Waves," *Journal of Computational Physics* 318 (2016): 1-21, https://doi.org/10.1016/j.jcp.2016.04.060.

5. U. Bosi, A. Engsig-Karup, C. Eskilsson, and M. Ricchiuto, "A Spectral/Hp Element Depth-Integrated Model for Nonlinear Wave-Body Interaction," *Computer Methods in Applied Mechanics and Engineering* 348 (2019): 222-249, https://doi.org/10.1016/j.cma.2019.01.020.

6. M. Göteman, M. Giassi, J. Engström, and J. Isberg, "Advances and Challenges in Wave Energy Park Optimization - A Review," *Frontiers in Energy Research* 8 (2020): 1-26, https://doi.org/10.3389/fenrg.2020.00026.

7. D. Golbaz, R. Asadi, E. Amini, et al., "Layout and Design Optimization of Ocean Wave Energy Converters: A Scoping Review of State-Of-The-Art Canonical, Hybrid, Cooperative, and Combinatorial Optimization Methods," *Energy Reports* 8 (2022): 15446-15479, https://doi.org/10.1016/j.egyr.2022.10.403.

8. T. Sayadi, P. J. Schmid, F. Richecoeur, and D. Durox, "Parametrized Data-Driven Decomposition for Bifurcation Analysis, With Application to Thermo-Acoustically Unstable Systems," *Physics of Fluids* 27, no. 3 (2015): 037102.

9. S. Roy, J. C. Hua, W. Barnhill, G. H. Gunaratne, and J. R. Gord, "Deco-convolution of Reacting-Flow Dynamics Using Proper Orthogonal and Dynamic Mode Decompositions," *Physical Review E* 91 (2015): 013001.

10. A. Treuille, A. Lewis, and Z. Popovic, "Model Reduction for Real-Time Fluids," *ACM Transactions on Graphics* 25 (2006): 826-834, https://doi.org/10.1145/1141911.1141962.

11. B. Liu, G. Mason, J. Hodgson, Y. Tong, and M. Desbrun, "Model-Reduced Variational Fluid Simulation," *ACM Transactions on Graphics* 34, no. 6 (2015): 1-12, https://doi.org/10.1145/2816795.2818130.

12. C. W. Rowley and S. T. Dawson, "Model Reduction for Flow Analysis and Control," *Annual Review of Fluid Mechanics* 49, no. 1 (2017): 387-417, https://doi.org/10.1146/annurev-fluid-010816-060042.

13. T. Lassila, A. Manzoni, A. Quarteroni, and G. Rozza, *Model Order Reduction in Fluid Dynamics: Challenges and Perspectives* (Cham: Springer International Publishing, 2014), 235-273.

14. S. S. Ravindran, "A Reduced-Order Approach for Optimal Control of Fluids Using Proper Orthogonal Decomposition," *International Journal for Numerical Methods in Fluids* 34, no. 5 (2000): 425-448.

15. A. Barbagallo, D. Sipp, and P. J. Schmid, "Closed-Loop Control of an Open Cavity Flow Using Reduced-Order Models," *Journal of Fluid Mechanics* 641 (2009): 1-50.

16. M. Strazzullo, F. Ballarin, R. Mosetti, and G. Rozza, "Model Reduction for Parametrized Optimal Control Problems in Environmental Marine Sciences and Engineering," *SIAM Journal on Scientific Computing* 40, no. 4 (2018): B1055-B1079, https://doi.org/10.1137/17M1150591.

17. A. Manzoni, A. Quarteroni, and G. Rozza, "Shape Optimization for Viscous Flows by Reduced Basis Methods and Free-Form Deformation," *International Journal for Numerical Methods in Fluids* 70, no. 5 (2012): 646-670.

18. A. Hay, J. Borggaard, I. Akhtar, and D. Pelletier, "Reduced-Order Models for Parameter Dependent Geometries Based on Shape Sensitivity Analysis," *Journal of Computational Physics* 229, no. 4 (2010): 1327-1352.

19. T. Bui-Thanh, K. Willcox, and O. Ghattas, "Parametric Reduced-Order Models for Probabilistic Analysis of Unsteady Aerodynamic Applications," *AIAA Journal* 46, no. 10 (2008): 2520-2529.

20. Z. Wang, I. Akhtar, J. Borggaard, and T. Iliescu, "Proper Orthogonal Decomposition Closure Models for Turbulent Flows: A Numerical Comparison," *Computer Methods in Applied Mechanics and Engineering* 237-240 (2012): 10-26.

21. J. Weller, E. Lombardi, M. Bergmann, and A. Iollo, "Numerical Methods for Low-Order Modeling of Fluid Flows Based on POD," *International Journal for Numerical Methods in Fluids* 63, no. 2 (2010): 249-268.

22. S. E. Ahmed, S. Pawar, O. San, A. Rasheed, T. Iliescu, and B. R. Noack, "On Closures for Reduced Order Models - A Spectrum of First-Principle to Machine-Learned Avenues," *Physics of Fluids* 33, no. 9 (2021): 091301.

23. T. Lassila, A. Manzoni, A. Quarteroni, and G. Rozza, "A Reduced Computational and Geometrical Framework for Inverse Problems in Hemodynamics," *International Journal for Numerical Methods in Biomedical Engineering* 29, no. 7 (2013): 741-776.

24. S. Sarmast, R. Dadfar, R. F. Mikkelsen, et al., "Mutual Inductance Instability of the Tip Vortices Behind a Wind Turbine," *Journal of Fluid Mechanics* 755 (2014): 705-731.

25. D. Xiao, Z. Lin, F. Fang, et al., "Non-intrusive Reduced Order Modeling for Multiphase Porous Media Flows Using Smolyak Sparse Grids," *International Journal for Numerical Methods in Fluids* 83, no. 2 (2016): 205-219, https://doi.org/10.1002/fld.4263.

26. K. Carlberg, C. Farhat, J. Cortial, and D. Amsallem, "The GNAT Method for Nonlinear Model Reduction: Effective Implementation and Application to Computational Fluid Dynamics and Turbulent Flows," *Journal of Computational Physics* 242 (2013): 623-647, https://doi.org/10.1016/j.jcp.2013.02.028.

27. K. Carlberg, M. Barone, and H. Antil, "Galerkin v. Least-Squares Petrov-Galerkin Projection in Nonlinear Model Reduction," *Journal of Computational Physics* 330 (2017): 693-734, https://doi.org/10.1016/j.jcp.2016.10.033.

28. I. Kalashnikova, M. Barone, and M. Brake, "A Stable Galerkin Reduced Order Model for Coupled Fluidstructure Interaction Problems," *International Journal for Numerical Methods in Engineering* 95, no. 2 (2013): 121-144, https://doi.org/10.1002/nme.4499.

29. B. Xu, H. Gao, M. Wei, and J. Hrynuk, "Global POD-Galerkin ROMs for Fluid Flows With Moving Solid Structures," *AIAA Journal* 60, no. 3 (2022): 1400-1414, https://doi.org/10.2514/1.J060795.

30. T. A. Brenner, R. L. Fontenot, P. G. Cizmas, and O. B. RW, "A Reduced-Order Model for Heat Transfer in Multiphase Flow and Practical Aspects of the Proper Orthogonal Decomposition," *Computers & Chemical Engineering* 43 (2012): 68-80, https://doi.org/10.1016/j.compchemeng.2012.04.003.

31. J. Li, T. Zhang, S. Sun, and B. Yu, "Numerical Investigation of the POD Reduced-Order Model for Fast Predictions of Two-Phase Flows in Porous Media," *International Journal of Numerical Methods for Heat & Fluid Flow* 29, no. 11 (2019): 4167-4204, https://doi.org/10.1108/HFF-02-2019-0129.

32. B. Battisti, T. Blickhan, G. Enchery, V. Ehrlacher, D. Lombardi, and O. Mula, "Wasserstein Model Reduction Approach for Parametrized Flow Problems in Porous Media," *ESAIM: Proceedings and Surveys* 73 (2023): 28-47, https://doi.org/10.1051/proc/202373028.

33. J. Barnett and C. Farhat, "Quadratic Approximation Manifold for Mitigating the Kolmogorov Barrier in Nonlinear Projection-Based Model Order Reduction," *Journal of Computational Physics* 464 (2022): 111348, https://doi.org/10.1016/j.jcp.2022.111348.

34. J. Barnett, C. Farhat, and Y. Maday, "Neural-Network-Augmented Projection-Based Model Order Reduction for Mitigating the Kolmogorov Barrier to Reducibility," *Journal of Computational Physics* 492 (2023): 112420, https://doi.org/10.1016/j.jcp.2023.112420.

35. C. Czech, M. Lesjak, C. Bach, and F. Duddeck, "Data-Driven Models for Crashworthiness Optimisation: Intrusive and Non-intrusive Model
Order Reduction Techniques,” *Structural and Multidisciplinary Optimization* 65, no. 190 (2022): 1–15, https://doi.org/10.1007/s00158-022-03282-1.

36. S. Aliyar, G. Ducrozet, B. Bouscasse, F. Bonnefoy, V. Sriram, and P. Ferrant, “Numerical Coupling Strategy Using HOS-OpenFOAM-MoorDyn for OC3 Hywind SPAR Type Platform,” *Ocean Engineering* 263 (2022): 112206, https://doi.org/10.1016/j.oceaneng.2022.112206.

37. T. Verbrugghe, V. Stratigaki, C. Altomare, J. M. Domínguez, P. Troch, and A. Kortenhaus, “Implementation of Open Boundaries Within a Two-Way Coupled SPH Model to Simulate Nonlinear Wavestructure Interactions,” *Energies* 12, no. 4 (2019), https://doi.org/10.3390/en12040697.

38. A. Mota, I. Tezaur, and G. Phlipot, “The Schwarz Alternating Method for Transient Solid Dynamics,” *International Journal for Numerical Methods in Engineering* 123, no. 21 (2022): 5036–5071, https://doi.org/10.1002/nme.6982.

39. I. Prusak, M. Nonino, D. Torlo, F. Ballarin, and G. Rozza, “An Optimisation-Based Domain-Decomposition Reduced Order Model for the Incompressible Navier-Stokes Equations,” *Computers & Mathematics With Applications* 151 (2023): 172–189, https://doi.org/10.1016/j.camwa.2023.09.039.

40. A. de Castro, P. Kuberry, I. Tezaur, and P. Bochev, “A novel partitioned approach for reduced order model–finite element model (ROM-FEM) and ROM-ROM coupling,” in *Earth and Space 2022* (2022): 475–489.

41. C. Hirt and B. Nichols, “Volume of Fluid (VOF) Method for the Dynamics of Free Boundaries,” *Journal of Computational Physics* 39, no. 1 (1981): 201–225, https://doi.org/10.1016/0021-9991(81)90145-5.

42. M. Bergmann, “Numerical Modeling of a Self-Propelled Dolphin Jump out of Water,” *Bioinspiration & Biomimetics* 17, no. 6 (2022): 065010, https://doi.org/10.1088/1748-3190/ac8fc8.

43. Siemens Digital Industries Software (Simcenter STAR-CCM+, version 2021.1.; Siemens, 2021).

44. G. Jin and M. Braza, “A Nonreflecting Outlet Boundary Condition for Incompressible Unsteady Navier-Stokes Calculations,” *Journal of Computational Physics* 107, no. 2 (1993): 239–253, https://doi.org/10.1006/jcph.1993.1140.

45. C. H. Bruneau and P. Fabrie, “Effective Downstream Boundary Conditions for Incompressible Navier-Stokes Equations,” *International Journal for Numerical Methods in Fluids* 19, no. 8 (1994): 693–705, https://doi.org/10.1002/fld.1650190805.

46. M. Bergmann, A. Ferrero, A. Iollo, E. Lombardi, A. Scardigli, and H. Telib, “A Zonal Galerkin-Free POD Model for Incompressible Flows,” *Journal of Computational Physics* 352 (2018): 301–325, https://doi.org/10.1016/j.jcp.2017.10.001.

47. S. Riffaud, M. Bergmann, C. Farhat, S. Grimberg, and A. Iollo, “The DGDD Method for Reduced-Order Modeling of Conservation Laws,” *Journal of Computational Physics* 437 (2021): 110336, https://doi.org/10.1016/j.jcp.2021.110336.

48. J. Borée, “Extended Proper Orthogonal Decomposition: A Tool to Analyse Correlated Events in Turbulent Flows,” *Experiments in Fluids* 35, no. 2 (2003): 188–192, https://doi.org/10.1007/s00348-003-0656-3.

49. B. R. Noack, K. Afanasiev, M. Morzynski, G. Tadmor, and F. Thiele, “A Hierarchy of Low-Dimensional Models for the Transient and Post-Transient Cylinder Wake,” *Journal of Fluid Mechanics* 497 (2003): 335–363, https://doi.org/10.1017/S0022112003006694.

50. M. Bergmann and L. Cordier, “Optimal Control of the Cylinder Wake in the Laminar Regime by Trust-Region Methods and POD Reduced-Order Models,” *Journal of Computational Physics* 227, no. 16 (2008): 7813–7840, https://doi.org/10.1016/j.jcp.2008.04.034.

51. J. L. Lumley, “The structure of inhomogeneous turbulence,” in *Atmospheric Turbulence and Wave Propagation*, eds. A. M. Yaglom and V. I. Tatarski (Moscow: Nauka, 1967), 166–178.

52. P. Holmes, J. L. Lumley, G. Berkooz, and C. W. Rowley, *Turbulence, Coherent Structures, Dynamical Systems and Symmetry* (Cambridge: Cambridge Monographs on Mechanics, 1996).

53. L. Cordier and M. Bergmann, *Proper Orthogonal Decomposition: An Overview*. Lecture series 2002-04, 2003-03 and 2008-01 on post-processing of experimental and numerical data, Von Karman Institute for Fluid Dynamics, 2008, VKI (Von Kármán Institute for Fluid Dynamics, 2008), 46.

54. L. Sirovich, “Turbulence and the Dynamics of Coherent Structures,” *Quarterly of Applied Mathematics* XLV, no. 3 (1987): 561–590.

55. B. Battisti, *Multi-Fidelity Multi-Scale Numerical Modelling of Wave Energy Converter Farms*. PhD thesis. (de Bordeaux: Politecnico di Torino & Université, 2024), https://theses.hal.science/tel-04633747.

56. M. Barrault, Y. Maday, N. Nguyen, and A. Patera, “An Empirical Interpolation Method: Application to Efficient Reduced-Basis Discretization of Partial Differential Equations,” *Comptes Rendus Mathematique* 339 (2004): 667–672, https://doi.org/10.1016/j.crma.2004.08.006.

57. C. Farhat, T. Chapman, and P. Avery, “Structure-Preserving, Stability, and Accuracy Properties of the Energy-Conserving Sampling and Weighting Method for the Hyper Reduction of Nonlinear Finite Element Dynamic Models,” *International Journal for Numerical Methods in Engineering* 102, no. 5 (2015): 1077–1110, https://doi.org/10.1002/nme.4820.

58. I. López, J. Andreu, S. Ceballos, I. Martínez de Alegría, and I. Kortabarria, “Review of Wave Energy Technologies and the Necessary Power-Equipment,” *Renewable and Sustainable Energy Reviews* 27 (2013): 413–434, https://doi.org/10.1016/j.rser.2013.07.009.