# High-order finite difference solution for 3D nonlinear wave-structure interaction

G. Ducrozet, H.B. Bingham, A.P. Engsig-Karup, P. Ferrant

► To cite this version:

G. Ducrozet, H.B. Bingham, A.P. Engsig-Karup, P. Ferrant. High-order finite difference solution for 3D nonlinear wave-structure interaction. Journal of Hydrodynamics, 2010, Proceedings of the 9th International Conference on Hydrodynamics, 22 (Suppl 1), pp.225-230. 10.1016/S1001-6058(09)60198-0. hal-01145092

HAL Id: hal-01145092

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

Submitted on 2 Jun 2025

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>

Distributed under a Creative Commons Attribution 4.0 International License
# High-order finite difference solution for 3D nonlinear wave-structure interaction

G. Ducrozet¹,², H.B. Bingham¹, A.P. Engsig-Karup³ & P. Ferrant²

¹ Department of Mechanical Engineering
Technical University of Denmark
Kgs Lyngby, Denmark
E-mail: hbb@mek.dtu.dk

² Laboratoire de Mécanique des Fluides
Ecole Centrale de Nantes
Nantes, France
E-mail: guillaume.ducrozet@ec-nantes.fr
pierre.ferrant@ec-nantes.fr

³ Department of Informatics and Mathematical Modeling
Technical University of Denmark
Kgs Lyngby, Denmark
E-mail: apek@imm.dtu.dk

## ABSTRACT

This contribution presents our recent progress on developing an efficient fully-nonlinear potential flow model for simulating 3D wave-wave and wave-structure interaction over arbitrary depths (i.e. in coastal and offshore environment). The model is based on a high-order finite difference scheme OceanWave3D presented in [1, 2]. A nonlinear decomposition of the solution into incident and scattered fields is used to increase the efficiency of the wave-structure interaction problem resolution. Application of the method to the diffraction of nonlinear waves around a fixed, bottom mounted circular cylinder are presented and compared to the fully nonlinear potential code XWAVE as well as to experiments.

**KEY WORDS:** High-Order Finite Differences; Nonlinear Decomposition; Scattering; OceanWave3D; XWAVE

## INTRODUCTION

Free surface flow problems in the naval and ocean engineering contexts are a subject of major interest. Particularly the need to accurately compute nonlinear wave loads acting on objects in the ocean environment is increasingly important. In the last decades there have been significant developments of several numerical methods addressing this problem ranging from potential flow solvers formulated in the frequency domain (with development to second or third order in nonlinearity) to fully nonlinear time-domain potential flow solvers up to Reynolds-Averaged Navier-Stokes Equations (RANSE) solvers. However, the latter

method usually involves computational efforts which impose severe limitations on domain size and/or resolution.

This study deals with the development of a fully-nonlinear potential flow model for 3D wave-wave and wave-structure interactions over arbitrary depths (i.e. able to solve the problem in coastal and offshore environment). The model is based on a high-order finite difference scheme presented in [1, 2]. The important properties of efficiency, accuracy and robustness for a model of this kind for wave propagation and wave/bottom interactions has been demonstrated. The efficient resolution of wave-structure interactions is the next step of interest in order to ensure that this computational tool is suitable for large-scale 3D nonlinear simulations. The XWAVE model [3] is used as a reference fully nonlinear potential model for which efficiency, accuracy and robustness have also been demonstrated. This model is based on a fully nonlinear Mixed Euler-Lagrange (MEL) formulation where a boundary element method (BEM) is used for the solution of the boundary integral equation formulation of the problem.

Firstly, the numerical methods used to solve the problem are briefly described. The general potential formulation is reviewed as well as the numerical approach adopted in the OceanWave3D model which is based on high-order finite differences. Following [3], we adopt a nonlinear decomposition of the wave-structure interaction solution into an incident and a scattered wave field. The incident part represents
the incoming wave which is known explicitly while the scattered part will take into account the presence of the structure inside the given incident wave field. Then, the MEL-BEM code XWAVE is also briefly presented. A first and preliminary comparison in terms of computational effort between the two models mentioned is presented. Finally, the Ocean-Wave3D model is applied to the nonlinear diffraction of regular long waves around a fixed, bottom-mounted circular cylinder. Comparisons with experiments and the XWAVE code are presented and discussed.

# NUMERICAL METHODS

## Formulation

We consider a cartesian coordinate system with the z axis defining the vertical and oriented upwards, and the xy plane (i.e. level $z = 0$) corresponds to the mean water level. The notation $\mathbf{x}$ stands for the $(x, y)$ vector. The fluid simulated is water, which is assumed to be incompressible and inviscid. The flow is also considered irrotational. With these assumptions one can apply potential flow theory, the velocity $\mathbf{V}$ derives from a potential $\mathbf{V}(\mathbf{x}, z, t) = \tilde{\nabla}\phi(\mathbf{x}, z, t)$, $\tilde{\nabla}$ representing the gradient and $\phi$ the scalar velocity potential function. Then, the continuity equation $(\tilde{\nabla} \cdot \mathbf{V} = 0)$ becomes the Laplace equation $\Delta\phi = 0$. On the free surface, the elevation is supposed to be single-valued and thus described by $z = \eta(\mathbf{x}, t)$. The slip condition and the continuity of the pressure give respectively the kinematic and dynamic free surface boundary conditions which can be used to account for the time-evolution of the free-surface boundary.

## OceanWave3D model

Following Zakharov [4], the fully-nonlinear free surface boundary conditions can be written in terms of surface quantities, namely $\eta$ and the surface potential $\tilde{\phi}(\mathbf{x}, t) = \phi(\mathbf{x}, \eta, t)$

$$ \frac{\partial \tilde{\phi}}{\partial t} = -g\eta - \frac{1}{2} |\nabla \tilde{\phi}|^2 + \frac{1}{2} (1 + |\nabla \eta|^2) \tilde{w}^2 \quad (1) $$

$$ \frac{\partial \eta}{\partial t} = (1 + |\nabla \eta|^2) \tilde{w} - \nabla \tilde{\phi} \cdot \nabla \eta \quad (2) $$

on $z = \eta(\mathbf{x}, t)$, with $\tilde{w}$ the vertical velocity at the free surface $\tilde{w}(\mathbf{x}, t) = \partial_z\phi(\mathbf{x}, \eta, t)$. To advance the unknowns $(\eta, \phi)$ in time, a classical explicit 4th order Runge-Kutta scheme can be used.

To account for the time evolution of the quantities of interest $(\eta, \tilde{\phi})$, one only need to evaluate this vertical

velocity at the free surface. This is done by solving the following Laplace problem for the velocity potential

$$ \phi = \tilde{\phi}, \quad z = \eta \quad (3) $$

$$ \nabla^2 \phi + \phi_{zz} = 0, \quad -h \le z < \eta \quad (4) $$

$$ (\mathbf{n}, n_z) \cdot (\nabla, \partial_z) \phi = 0, \quad (\mathbf{x}, z) \in \partial\Omega \quad (5) $$

with (5) the kinematic boundary conditions of the different boundaries $\partial\Omega$ of the finite domain, including the bottom whose profile can be chosen arbitrarily.

For the direct solution of this problem, a $\sigma$-transformation is employed in the vertical. At the same time, the description of fixed structures extending in the vertical from the bottom is made possible by the use of curvilinear boundary-fitted meshes. The physical domain can then be mapped to a simple computational domain, which leads to an efficient solution strategy. The bathymetry can be chosen freely as well as the geometry of the vertical bottom-mounted structure studied in the domain. This resulting model is suitable for the simulation of wave-wave, wave-bottom and wave-structure interactions.

The spatial derivatives operators are approximated by central high-order finite differences approximations except at the boundaries where off-centered stencils are employed. The order of the scheme is left arbitrary while the mesh is defined as logically structured and non-uniform. The resulting coupled system of equations is sparse and unsymmetric, thus we solve using the classical GMRES iterative method. Here, the preconditioning strategy is crucial to obtaining an efficient solution (a small number of iterations). As demonstrated in [2], the use of multigrid preconditioning makes it possible to achieve an optimal scaling of the overall solution effort (i.e. linear scaling with the total number of grid points).

More details about the OceanWave3D method can be found in Bingham and Zhang [1] for the 2D case and Engsig-Karup et al. [2] for the 3D case. In particular, stability and accuracy analyses in both two and three space dimensions have been performed together with a range of validation tests that confirm that the model is robust, efficient and high-order accurate.

## Nonlinear decomposition of the solution

In order to treat efficiently the wave-structure interaction problem, a fully nonlinear decomposition is
used. The whole problem is split up into an inci- dent part describing the undisturbed wave field and a scattered component taking care of the influence of the structure inside the domain. This approach has been widely used and validated (introduced first [3] in the XWAVE potential model and [5, 6] in the SWENSE solver based on RANS Equations). The solution of the whole problem ($\eta, \phi$) gives

$$ \eta = \eta^I + \eta^S, \quad \phi = \phi^I + \phi^S \quad (6) $$

with ($\eta^S, \phi^S$) the scattered components and ($\eta^I, \phi^I$) the incident wavefield which is known explicitly from a dedicated wave-solver (e.g. here we use stream function theory [7]). The fully nonlinear decomposition (6) is then inserted into the free surface boundary conditions (1) and (2) as well as other boundary conditions. These are reformulated in terms of scattered quantities (i.e. evolution equations $\partial_t\eta^S$ and $\partial_t\phi^S$). The incident components will then act as forcing terms into these modified boundary conditions. Note that since the problem is fully nonlinear, the modified equations (1) and (2) have to be satisfied at the exact free surface position and thus, the incident wave field will often be evaluated above its surface of definition. However, typical potential flow based wave theories (such as stream function theory) allow for an analytic continuation of the incident potential above the true wave elevation.

The advantages of this procedure are: i) the incident wave is treated explicitly by a dedicated model and is thus not altered during its propagation from the outer boundaries up to the body of interest; ii) the scattered waves are mainly radial which allows the use of a refined grid close to the body of interest and coarse grids far from it. This will result in a considerable reduction of computational effort during the resolution of the wave-structure interaction problem.

## XWAVE model

For comparison with the OceanWave code, the XWAVE model is used [3]. This model has been widely validated in the context of fully nonlinear time domain diffraction and radiation test cases (e.g see [3, 8]). This is based on a Semi-Lagrangian formalism in which the mesh on the free surface is only left free in the vertical direction. A classical Boundary Integral Equation Method (BIEM) is used to transform the 3D problem onto a 2D one formulated on the boundary of the domain (through the use of Green's theorem). A Boundary Element Method (BEM) is then used

for the solution of this BIEM problem. The resolution scheme leads to a full, non-symmetric square system of equations to be solved for the singularity distributions at the boundaries of the computational domain. This is done by classical iterative GMRES procedure with diagonal preconditioning. Notice that the previously introduced nonlinear decomposition is also used in the XWAVE model in which it has shown its applicability and efficiency.

## Computational effort

In the fully nonlinear model XWAVE, the boundary geometry evolves at each time step requiring both the construction and the solution of a new linear system at each iteration in time. Preconditioned GMRES is used to solve this system requiring a computational effort evolving as $O(N^2)$. At the same time, it has been demonstrated in [2], that for OceanWave3D, the use of multigrid preconditioning makes it possible to achieve a linear scaling of the overall solution effort with the number of points (i.e. solution obtained in $O(N)$).

In order to compare both approaches, one has to keep in mind that the BIEM uses Green's theorem to transform the 3D problem onto a 2D problem on the domain boundaries. The number of points in both methods is there not equivalent, OceanWave3D discretizing the whole fluid volume. However, the efficiency of the high-order finite difference-based model allow the use of very few grid points in the vertical to achieve good accuracy (typically at most 10 nodes clustered in the vertical toward the free surface) and thereby makes it competitive. The additional expense of the volume discretization is thus not so problematic compared to a surface one (such as the XWAVE code).
be faster.

However, note that this quick comparison is just to give an idea of the respective computational effort of both methods. More in-depth comparisons have still to be performed in order to compare at given accuracy, the effective computational effort of each of the method. Furthermore, for now the meshes used are not strictly equivalent between both models which prevents an exact comparison.

## NONLINEAR DIFFRACTION

In this part, nonlinear diffraction by a cylinder in long waves will be studied in a domain of small finite depth. In particular, the wave run-up around the cylinder is compared to the experiments conducted by Kriebel [9]. In order to assess the accuracy of the OceanWave3D model, the results will also be compared to the fully nonlinear potential model XWAVE.

### Problem set-up

<figure><img src="image_2.png" /><figcaption>Figure 2: Scheme of the computational domain</figcaption></figure>

Three of the cases given in the paper by Kriebel [9] are reproduced. The interaction of long regular waves (wavelength $\lambda$, wavenumber $k = 2\pi/\lambda$) of different amplitudes are studied. The domain is at constant depth $h$ with $kh = 1.036$ (i.e. rather small depth) and the cylinder radius $R$ is small compared to the wavelength $kR = 0.374$. The influence of nonlinearities is studied with the three cases chosen in a range where no wave-breaking occurs during the experiments: the wave steepnesses of incident wave field (of amplitude $A$) are $kA = 0.061$, $kA = 0.1025$ and $kA = 0.143$. Figure 2 presents an example of the domain of resolution used in OceanWave3D (similar to the one used in XWAVE). Only half of the domain is necessary to solve the whole problem and notice that on this sketch the vertical direction is amplified by a factor 3 for clarity (this to remind that the domain is quite

<figure><img src="image_3.png" /><figcaption>Figure 1: Scaling of computational effort for XWAVE and OceanWave3D with number of points $N$ on the surface of the computational domain</figcaption></figure>

An example of the scaling of the solution effort for both the OceanWave3D and the XWAVE methods is provided on Fig. 1 where the total CPU time to compute 100 time steps is presented. The scaling is done with respect to the number of points $N$ at the surface of the domain in order to make a fair comparison between the volume code and the surface one. Results are obtained for a nonlinear diffraction test-case such as the ones presented in next section. Typical parameters are used in terms of GMRES tolerance or time-stepping in both methods. A 6th order finite difference scheme is used in OceanWave3D.

Figure 1 presents the two reference curves $O(N^2)$ and $O(N')$. Notice that the $O(N')$ curve is obtained with the total number of points in the volumic domain $N'$ (i.e. not directly from the number of points $N$). Indeed, for OceanWave3D, the linear scaling of the computational effort is obtained with respect to the total number of points in the domain.

Once this is noted, one can see that both of models behave like expected: XWAVE is $O(N^2)$ and OceanWave3D $O(N')$. The interesting feature is that for a given number of points at the surface of the computational domain, the high-order finite difference scheme is more efficient than the MEL-BEM one. The number of points is here typical for a nonlinear diffraction case. The XWAVE code is expected to be more efficient for a very few number of points (taking into account the different behaviours $O(N^2)$ vs. $O(N')$) while for typical situation, OceanWave3D will
shallow).

Concerning the numerical resolution of the problem, in both models it has been chosen to solve precisely up to the 3rd harmonic. This imposes a number of points per 3rd harmonic wavelength of around 16 to ensure a precise resolution, i.e. around 100 points per $λ$. The OceanWave3D code uses here 8 points in the vertical and a 6th order finite difference scheme. The time integration is done with a classical 4th order Runge-Kutta scheme in both models with a Courant number $C_r = 0.6$.

## Results

On the following figures are compared the fully nonlinear simulations (OceanWave3D and XWAVE) with Stokes 2nd order theory and experiments made by Kriebel [9]. Figures 3, 4 & 5 stand for $kA = 0.061$, $kA = 0.1025$ & $kA = 0.143$ respectively.

<figure><img src="image_4.png" /><figcaption>Figure 3: Run-up along cylinder, $kh = 1.036$, $kR = 0.374$, $kA = 0.061$</figcaption></figure>

The first test case presented on Fig. 3 deals with the interaction of low steepness regular waves with the vertical cylinder ($kA = 0.061$). Both fully nonlinear simulations are in excellent concordance and in fairly good agreement with experiments. The results of the 2nd order theory are satisfactory with some discrepancies observed in the back of the cylinder ($\theta > 2\pi/3$). It is then noticeable that even with low steepness waves, high order nonlinear effects (superior than 2) occur during the diffraction around this cylinder. These effects appear to be not negligible

and correctly captured by OceanWave3D and XWAVE.

The amplitude of the incoming waves is then enlarged to see the evolution of nonlinear phenomena occuring in the diffraction process. Figure 4 presents the result of the run-up around the cylinder for an incoming wave of steepness $kA = 0.1025$. Clear differences are now observed between Stokes 2nd order theory and experiments. High-order effects are then essential in the modelization of such phenomenon. The accuracy of the comparisons between both fully nonlinear methods and experiments shows their ability to deal with highly nonlinear free surface flows. Notice that slight filtering is applied in both models to take care of small instabilities which may appear close to the cylinder, particularly to the rear.

<figure><img src="image_5.png" /><figcaption>Figure 4: Run-up along cylinder, $kh = 1.036$, $kR = 0.374$, $kA = 0.1025$</figcaption></figure>

Finally, Fig. 5 presents the results for an incoming wave steepness $kA = 0.143$. Previous conclusions are confirmed with a 2nd order theory which greatly underpredicts the wave run-up in the front of the cylinder. At the same time, OceanWave3D captures precisely most of the physics involved in this diffraction case (and particularly the high-order effects which are not part of the 2nd order theory). The XWAVE code is still in good agreement with experiments with slight discrepancies observed (small attenuation of the run up at the front of the cylinder for instance). These are probably due to the smoothing process necessary to the stability of the solution.
<figure><img src="image_6.png" /><figcaption>Figure 5: Run-up along cylinder, $kh = 1.036$, $kR = 0.374$, $kA = 0.143$</figcaption></figure>

## CONCLUSION

The high-order finite differences numerical model *Ocean Wave 3D* has been presented. An approach based on a nonlinear decomposition of incident and scattered wave field is used to enhance the efficiency of the wave-structure interaction problem resolution. Preliminary comparisons in terms of computational effort between *Ocean Wave 3D* and the fully nonlinear MEL-BEM model *XWAVE* are exposed. The high-order finite differences model appears to be more efficient for typical situation and thus suitable for large-scale simulations of wave-wave, wave-bottom and wave-structure interactions. Finally, a validation on the case of the nonlinear wave diffraction of a vertical bottom-mounted cylinder in long regular waves is detailed. Precise comparisons with experiments and the *XWAVE* model are provided in a wide range of nonlinearities which point out the accuracy and robustness of our numerical tool.

The inclusion of fixed and floating bodies with complex geometries is a topic of ongoing development. Future work also includes the extension to irregular wave fields through the use of High-Order Spectral method (see e.g. Ducrozet et al. [10]). This will allow the description of any complex incident wave fields and their interactions with structures in domains with

## ACKNOWLEDGMENTS

The authors wish to thank Délégation Générale de l'Armement (project PDE 08C0039) for partially funding this work, and the Danish Center for Scientific Computing for supercomputing resources.

## References

[1] H.B. Bingham and H. Zhang. On the accuracy of finite-difference solutions for nonlinear water waves. *J. Eng. Math.*, 58:211–228, 2007.

[2] A.P. Engsig-Karup, H.B. Bingham, and O. Lindberg. An efficient flexible-order model for 3d nonlinear water waves. *J. Comp. Phys.*, 228:2100–2118, 2009.

[3] P. Ferrant. Fully non-linear interactions of long-crested wave packets with a three-dimensional body. In *Proc. 22nd ONR Symposium on Naval Hydrodynamics*, pages 403–415, Washington, USA, 1998.

[4] V. Zakharov. Stability of periodic waves of finite amplitude on the surface of a deep fluid. *J. Appl. Mech. Tech. Phys.*, pages 190–194, 1968.

[5] P. Ferrant, L. Gentaz, B. Alessandrini, and D. Le Touzé. A potential/ranse approach for regular water wave diffraction about 2d structures. *Ship Technology Research*, 50(4), 2003.

[6] L. Gentaz, R. Luquet, B. Alessandrini, and P. Ferrant. Numerical simulation of the 3d viscous flow around a vertical cylinder in non-linear waves using an explicit incident wave model. In *Proc. 23rd Int. Conf. on Offshore Mech. and Arctic Engng.*, 2004.

[7] J.D. Fenton. The numerical solution of the steady water wave problem. *Comp. & Geosc.*, 14(3):357–368, 1988.

[8] B. Büchmann, P. Ferrant, and J. Skourup. Run-up on a body in waves and current. fully nonlinear and finite-order calculations. *Applied Ocean Research*, 22(6):349–360, 2000.

[9] D. L. Kriebel. Nonlinear wave interaction with a vertical circular cylinder. Part II: wave run-up. *Ocean Engng.*, 19(1):75–99, 1992.

[10] G. Ducrozet, F. Bonnefoy, D. Le Touze, and P. Ferrant. Implementation and validation of nonlinear wavemaker models in a HOS numerical wave tank. *International Journal of Offshore and Polar Engineering*, 16(3):161–167, SEP 2006.