<figure><img src="image_1.png" /></figure>

DTU Library

# A high-order shifted boundary method for water waves and floating bodies

Visbech, Jens; Engsig-Karup, Allan P.; Bingham, Harry B.; Amini-Afshar, Mostafa; Ricchiuto, Mario

**Published in:**
Proceedings of 39th International Workshop on Water Waves and Floating Bodies

**Publication date:**
2024

**Document Version**
Publisher's PDF, also known as Version of record

Link back to DTU Orbit

**Citation (APA):**
Visbech, J., Engsig-Karup, A. P., Bingham, H. B., Amini-Afshar, M., & Ricchiuto, M. (2024). A high-order shifted boundary method for water waves and floating bodies. In *Proceedings of 39th International Workshop on Water Waves and Floating Bodies IWWWFB*.

## General rights

Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright owners and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights.

* Users may download and print one copy of any publication from the public portal for the purpose of private study or research.

* You may not further distribute the material or use it for any profit-making activity or commercial gain

* You may freely distribute the URL identifying the publication in the public portal

If you believe that this document breaches copyright please contact us providing details, and we will remove access to the work immediately and investigate your claim.
# A high-order shifted boundary method for water waves and floating bodies

Jens Visbech¹, Allan P. Engsig-Karup¹, Harry B. Bingham², Mostafa Amini-Afshar², and Mario Ricchiuto³

¹ Dept. of Applied Mathematics, Technical University of Denmark, DK

² Dept. of Civil & Mechanical Eng., Technical University of Denmark, DK

³ Team CARDAMOM, INRIA, U. Bordeaux, CNRS, Bordeaux, France

Email to the corresponding author: jvis@dtu.dk

## INTRODUCTION

We consider the setting of fully nonlinear potential flow when simulating free surface waves using finite volume/cell/element methods. Here, the fluid domain is bounded by arbitrarily complex boundaries in topology and shape but also time-dependent (deforming or/and moving). Therefore, modeling such problems can be challenging and time-consuming from a meshing-point-of-view. If mesh updates are required [1], the domain should be re-meshed at each discrete time step, which provides a heavy computational burden. One way to circumvent this challenge is to consider unfitted/embedded/immersed boundary-type methods, where a simple – easy-to-generate – regular-shaped mesh is constructed once. Hereafter, the true domain is placed – and freely moved – on top of this background mesh. Enforcing boundary conditions is how unfitted methods differ from each other. One – relatively new – unfitted boundary method is the shifted boundary method (SBM), which we define more thoroughly in the upcoming section. Unfitted approaches have been developed for wave problems in recent years, where we highlight a few: The harmonic polynomial cell method with immersed type boundaries and/or overlapping grids in [2]. The same techniques have also been applied to the finite difference method in [3]. This abstract presents our recent developments toward modeling water waves and floating bodies using the SBM combined with a high-order spectral element method (SEM). Ultimately, we aim to model highly nonlinear wave-wave and wave-structure interactions using these techniques.

## SHIFTED BOUNDARY METHOD

Main & Scovazzi originally proposed the shifted boundary idea in [4, 5] for solving boundary value problems using the finite element method. The key idea behind the SBM is to embed the true domain, $\Omega \in \mathbb{R}^d$, where $d = \{1, 2, 3\}$, onto a regular-shaped affine mesh, denoted as the surrogate domain, $\bar{\Omega}_h \in \mathbb{R}^d$. Now, the problem is solved on $\bar{\Omega}_h$ and its surrogate boundary, $\bar{\Gamma}_h = \partial\bar{\Omega}_h$; however, the

boundary information is only defined on the true boundary, $\Gamma = \partial\Omega$. To close this problem, two ideas are combined: i) a mapping between points on the true boundary, $x = (x, y) \in \Gamma$ and points on the surrogate boundary, $\bar{x} = (\bar{x}, \bar{y}) \in \bar{\Gamma}_h$, as defined by

<figure><img src="image_2.png" /><figcaption>Figure 1: Concept of the SBM with $\Omega \subseteq \bar{\Omega}_h$.</figcaption></figure>

$$ \mathcal{M}(\bar{x}) : \bar{\Gamma}_h \mapsto \Gamma, \quad \bar{x} \mapsto x, \qquad (1) $$

and ii) a Taylor expansion to satisfy the boundary condition and conserve the optimal convergence rate of the chosen numerical method. The SBM naturally incorporates curved geometrical features, requires no re-meshing, and avoids the small-cut-cell problem that classical embedded cut-type methods suffer from. See Figure 1 for the SBM concept. Note that on the figure,
$\Omega \subseteq \bar{\Omega}_h$, yet $\bar{\Omega}_h \subseteq \Omega$, is also a possibility. Since the original work, many developments have been made for this numerical technique, including a high-order extension in [6]. The SBM and SEM connection was recently established for the two-dimensional Poisson problem in [7]. Here, a polynomial correction was employed to avoid the explicit evaluation of the Taylor series, which becomes cumbersome for high-order methods and three-dimensional problems [8].

## GOVERNING EQUATIONS

We consider the true fluid domain, $\Omega$, that is bounded by the time-depended free surface, $\Gamma^\eta$, and the bathymetry, $\Gamma^b$. Moreover, we have rigid walls, $\Gamma^w$, or periodic boundaries, $\Gamma^p$, both defined vertically at the ends of $\Omega$. Lastly, body boundaries are denoted by $\Gamma^{\text{body}}$. We adopt the setting of potential flow, such that the velocities, $u = (u, w)^T$, are expressed as $u = \nabla\phi$ where $\nabla = (\partial_x, \partial_y)^T$. Now, we seek to find $\phi \in C^2(\Omega)$ by solving the Laplace problem given as

$$ 
\begin{array}{ll}
\nabla^2 \phi = 0, & \text{in } \Omega \times \mathcal{T}, \\
\phi = \phi_\eta, & \text{on } \Gamma^\eta \times \mathcal{T}, \\
\nabla \phi \cdot \mathbf{n} = q_i, & \text{on } \Gamma^i, \\
\phi|_{x=x_{\text{min}}} = \phi|_{x=x_{\text{max}}}, & \text{on } \Gamma^p \text{ if } \Gamma^p \notin \emptyset,
\end{array}
\quad (2) $$

Figure 2: Embedding example of $\Gamma^\eta$, $\Gamma^b$, and $\Gamma^{\text{body}}$ in a background mesh and their respective normals, $\mathbf{n}$.

$$ \partial_t \eta = -\partial_x \eta \partial_x \phi_\eta + w_\eta (1 + \partial_x \eta \partial_x \eta), \quad \text{and} \quad \partial_t \phi_\eta = -g\eta - \frac{1}{2} \partial_x \phi_\eta \partial_x \phi_\eta + \frac{1}{2} w_\eta^2 (1 + \partial_x \eta \partial_x \eta). \quad (3) $$

<figure><img src="image_3.png" /></figure>

## WEAK FORMULATIONS

The boundary data in (2) are enforced weakly through an Aubin-type penalty formulation; see [7] and note that other formulations also can be used. For the conformal problem with $\Gamma^p \in \emptyset$, we get

$$ (\nabla\phi, \nabla v)_{\Omega} + \tau(\phi, v)_{\Gamma^\eta} - (\nabla\phi \cdot \mathbf{n}, v)_{\Gamma^\eta} = \tau(\phi_\eta, v)_{\Gamma^\eta} + (q_i, v)_{\Gamma^i}, \quad (4) $$

where $v$ is a test function and $(a, b)_S = \int_S a \cdot b \, dS$ is the inner product. Moreover, an implicit summation over $i \in \{\text{b, w, body}\}$ is used. Also, $\tau$ is a problem-depended penalty parameter. Now, we make the shift from $\Omega$ to $\bar{\Omega}_h$ and similarly for the boundaries. Therefore, we get a problem entirely defined in the surrogate domain as

$$ (\nabla\phi, \nabla v)_{\bar{\Omega}_h} + \tau(\phi, v)_{\bar{\Gamma}^\eta} - (\nabla\phi \cdot \bar{\mathbf{n}}, v)_{\bar{\Gamma}^\eta} = \tau(\bar{\phi}_\eta, v)_{\bar{\Gamma}^\eta} + (\bar{q}_i, v)_{\bar{\Gamma}^i}, \quad (5) $$

where $\bar{\mathbf{n}}$ is the outward facing normal on $\bar{\Gamma}_h$. The challenge is now to determine the value of the boundary data, $\bar{\phi}_\eta$ and $\bar{q}_i$, on $\bar{\Gamma}_h$ using $\phi_\eta$ and $q_i$ from $\Gamma$, while preserving consistency between the two. For this, we apply a polynomial correction (omitted for conciseness, see [7, 8]), as

$$ 
\begin{aligned}
& (\nabla\phi, \nabla v)_{\bar{\Omega}_h} + \tau(\phi(\mathbf{x}), v)_{\bar{\Gamma}^\eta} - (\nabla\phi \cdot \bar{\mathbf{n}}, v)_{\bar{\Gamma}^\eta} + \bar{n}_n(\nabla(\phi(\mathbf{x}) - \phi) \cdot \mathbf{n}, v)_{\bar{\Gamma}^i} - \bar{n}_t(\nabla\phi \cdot \mathbf{t}, v)_{\bar{\Gamma}^i} \\
& = \tau(\phi_\eta, v)_{\bar{\Gamma}^\eta} + \bar{n}_n(q_i, v)_{\bar{\Gamma}^i},
\end{aligned}
\quad (6) $$

where $\bar{n}_n = \bar{\mathbf{n}} \cdot \mathbf{n}$ and $\bar{n}_t = \bar{\mathbf{n}} \cdot \mathbf{t}$ with $\mathbf{t}$ being the tangent on $\Gamma$. Moreover, recall that $\mathbf{x} = \mathcal{M}(\bar{\mathbf{x}})$ which highlights the need for evaluating the basis – that is defined on $\bar{\Omega}_h$ – on $\Gamma$. Ultimately, this formulation allows us to simulate nonlinear free surface waves in a domain with an arbitrary curved/complex bathymetry and bodies with a deforming/moving free surface without re-meshing.
## NUMERICAL DISCRETIZATION

We adopt the usual *method of lines* approach, where the surrogate domain, $\bar{\Omega}_h$, is discretized using the SEM, see, e.g., [9] for the discrete approach to solve (2)-(3). We employ structured quadrilateral elements to represent $\bar{\Omega}_h$. The mapping, $x = M(\bar{x})$, is chosen to be purely vertical for $\Gamma^\eta$ and $\Gamma^b$ and along $n$ for $\Gamma^{\text{body}}$. See Figure 2 for a visualization of how $\Gamma^\eta$, $\Gamma^b$, and $\Gamma^{\text{body}}$ are embedded.

## RESULTS

**Verification of the Laplacian:** To verify the embedded Laplacian, we perform *p*-convergence studies of nonlinear stream functions waves with different wave numbers, $kh$, and percentages of maximum steepness, $\epsilon/\epsilon_{\text{max}}$, on a periodic domain. Here, $\epsilon = H/L$ is the wave steepness and $\epsilon_{\text{max}}$ is the maximum possible wave steepness for a given $kh$. With this, $\Gamma^w = \Gamma^{\text{body}} = \emptyset$ and $h$ is constant. In Figure 3, a study of a $kh = 4$ and $\epsilon/\epsilon_{\text{max}} = 70\%$ wave on three different meshes is seen, where $N_x$ and $N_z$ is the number of elements in the two directions. For this test, $\Gamma^b$ is embedded inside the bottom row of elements in $\bar{\Omega}_h$, whereas $\Gamma^\eta$ is located just above the top row. This results in both interpolation ($\Gamma^b \subseteq \bar{\Omega}_h$) and extrapolation ($\Gamma^\eta \not\subseteq \bar{\Omega}_h$) of the basis function as discussed in [7]. From the figure, exponential decay in the error is confirmed when increasing the polynomial order of the basis, although the curved free surface is represented on 2, 4, or 6 linear elements.

**Propagation of a stream function wave:** A periodic wave with $kh = 2\pi$ and $\epsilon/\epsilon_{\text{max}} = 70\%$ is propagated for 10 periods. The placement of $\Gamma^\eta$ and $\Gamma^b$ is as given above. At the final time, the free surface elevation is plotted against the true solution in Figure 4. The simulation is performed on a mesh with $(N_x, N_z) = (10, 2)$ for $P = 6$, and 1% of the energy in the top mode is filtered out between each time step. From the figure, the wave shape remains stable despite the embedded computation of the Laplacian. Minor numerical dispersive and diffusive errors are seen as expected for a high-order accurate numerical scheme for the resolution chosen.

<figure><img src="image_4.png" /><figcaption>Figure 3: Convergence study in $P$ for embedded Laplacian of periodic stream function wave with $kh = 4$ and $\epsilon/\epsilon_{\text{max}} = 70\%$ on three meshes.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: Stream function solution ($kh = 2\pi$ and $\epsilon/\epsilon_{\text{max}} = 70\%$) after 10 periods.</figcaption></figure>

**Horizontally forced cylinder in an infinite domain:** As a last numerical experiment, we displace a cylinder of radius $R = 1$ [m] horizontally by $U/\omega \sin(\omega t)$ in an infinite fluid. Here, $U = 0.2$ [m/s] and $\omega = 1$ [$s^{-1}$]. This gives an inhomogeneous Neumann condition on $\Gamma^{\text{body}}$ as $q_{\text{body}} = U \cos(\omega t)$. Dirichlet conditions based on potential flow theory are enforced on the remaining boundaries. A background mesh with elements of size $0.25 \times 0.25$ [m] and order $P = 4$ is constructed once. From this, the cylinder is moved horizontally due to the sinusoidal displacement, and the surrogate domain, $\bar{\Omega}_h$, is formed accordingly. About 24 linear surrogate boundary elements represent the cylinder depending on its location. The horizontal force signal, $F_x$, is obtained by integrating the linearized Bernoulli's equation over the body and validated against the true analytical solution as seen in Figure 5. The figure shows acceptable visual agreement despite the curved nature of the cylinder. With closer inspection (not visually
possible), small jumps in the force signal can be observed. This numerical artifact is generated by varying – mesh-dependent – levels of truncation errors, ultimately leading to discontinuities in the time signal. This phenomenon can also be seen in [3] and is subject to further investigation.

## CONCLUSION AND PERSPECTIVES

We have presented preliminary results for modeling water waves and bodies using a high-order shifted boundary method. Waves and forced body motions were modeled on very simple meshes, which still included curvature and circumvented re-meshing of the domain. These results indicate that this novel high-order embedded technique can be useful for modeling nonlinear wave-wave and wave-structure interactions. We aim to show more elaborate numerical cases for nonlinear wave-structure interaction at the workshop.

<figure><img src="image_6.png" /><figcaption>Figure 5: Force signal of a horizontally forced cylinder in an infinite domain.</figcaption></figure>

## ACKNOWLEDGMENTS

The presented research is carried out in connection with JV's Ph.D. project: "New advanced simulation techniques for wave energy converters". Funded as a DTU scholarship at DTU Compute. Part of the work was carried out during a visit to INRIA.

## REFERENCES

[1] Engsig-Karup, A. P., Monteserin, C., and Eskilsson, C. 2019. *A Mixed Eulerian–Lagrangian Spectral Element Method for Nonlinear Wave Interaction with Fixed Structures*. Water Waves 1(2), 315–342.

[2] Hanssen, F.-C. W., and Greco, M. 2021. *A potential flow method combining immersed boundaries and overlapping grids: Formulation, validation and verification*. Ocean Eng. 227, 108841.

[3] Xu, Y., Bingham, H. B., and Shao, Y. 2023. *A high-order finite difference method with immersed boundary treatment for fully nonlinear wave structure interaction*. Appl. Ocean Res. 134, 103535.

[4] Main, A., and Scovazzi, G. 2018. *The shifted boundary method for embedded domain computations. Part I: Poisson and Stokes problems*. J. Comput. Phys. 372, 972–995.

[5] Main, A., and Scovazzi, G. 2018. *The shifted boundary method for embedded domain computations. Part II: Linear advection diffusion and incompressible Navier Stokes equations*. J. Comput. Phys. 372, 996–1026.

[6] Atallah, N. M., Canuto, C., and Scovazzi, G. 2022. *The high-order Shifted Boundary Method and its analysis*. Comput. Methods Appl. Mech. Eng. 394, 114885.

[7] Visbech, J., Engsig-Karup, A. P., and Ricchiuto, M. 2023. *A spectral element solution of the Poisson equation with shifted boundary polynomial corrections: influence of the surrogate to true boundary mapping and an asymptotically preserving Robin formulation*. Preprint on arXiv (2310.17621).

[8] Ciallella, M., Gaburro, E., Lorini, M., and Ricchiuto, M. 2023. *Shifted boundary polynomial corrections for compressible flows: high order on curved domains using linear meshes*. Comput. Appl. Math. 441, 127698.

[9] Engsig-Karup, A. P., Eskilsson, C., and Bigoni, D. 2016. *A stabilised nodal spectral element method for fully nonlinear water waves*. J. Comput. Phys. 318, 1–21.