PaperlayerView exact revision on arXiv

Source-linked AI summary

A high-order polynomial-corrected shifted boundary method for simulating fully nonlinear water waves

Jens Visbech, Allan P. Engsig-Karup, Harry B. Bingham, Mario Ricchiuto

arXiv:2608.28123v1physics.flu-dynmath.NA

TL;DR

The paper addresses accurate simulation of fully nonlinear water waves with moving, deforming, and curved boundaries without remeshing. It introduces a high-order unfitted framework combining shifted-boundary spectral elements, polynomial-preserving gradient recovery, and arbitrary-order free-surface discretization, with verification and validation across nonlinear wave cases. The framework supports geometrically flexible simulations, while high-amplitude solitary-wave interaction remains bounded by free-surface topology changes.

  • Problem

    Unfitted discretizations for fully nonlinear potential-flow wave propagation and wave-structure interaction are relatively underexplored despite highly time-dependent curved and complex boundaries.

  • Method

    The framework combines polynomial-corrected shifted-boundary spectral elements with polynomial-preserving gradient recovery and arbitrary-order finite differences with hyperviscosity for the free surface.

  • Results

    The proposed framework provides high-order unfitted simulation of highly nonlinear water waves with moving and deforming geometry, supported by convergence, speed, long-time, bathymetry, and wall-interaction studies.

  • Takeaways & Limitations

    The model offers geometric flexibility for simulating nonlinear waves and their interactions with bathymetry and offshore structures on unfitted meshes.

  • Takeaways & Limitations

    For the most nonlinear high-amplitude solitary-wave case, the simulation breaks down when the free-surface elevation becomes non-single-valued.

Abstract

from arXiv · show

We present a novel unfitted computational framework for simulating fully nonlinear potential flow-based water waves. Focusing on wave propagation, we describe the core methodology, which involves a high-order polynomial-corrected shifted-boundary approximation on unfitted spectral elements. This approach allows for the simulation of a curved, highly time-dependent (moving and deforming) free surface affected by bathymetric changes. All that on a simple Cartesian background mesh without re-meshing or further approximations of non-affine geometric features. In addition, we highlight the importance of proper gradient recovery using a polynomial-preserving technique to accurately capture the vertical free-surface velocity. Ultimately, the goal is to develop a high-order convergent numerical scheme capable of simulating highly nonlinear waves over long periods. This is achieved through an arbitrary-order finite-difference approximation of the free surface combined with added hyperviscosity for numerical stability. We present verification and validation test cases for wave propagation in both periodic and finite domains. Emphasis is placed on convergence studies, the justification for using high-order approximations, the importance of optimal gradient recovery, long-time simulation of highly nonlinear waves, and nonlinear wave interactions with both affine and curved bathymetry changes, as well as with vertical walls, for highly nonlinear stream function waves and high-amplitude solitons. The novel computational model offers a comprehensive, all-in-one framework for simulating ocean waves and their interactions with offshore structures. It provides significant geometric flexibility, enabling boundaries such as the free surface to move and deform over time without the need to re-mesh a boundary-fitted mesh.

1 Introduction

The paper motivates high-order numerical models for nonlinear water waves because physical and numerical approximations constrain fidelity and long-time accuracy. It introduces an unfitted shifted-boundary framework for time-dependent free surfaces and reports verification and validation studies.

  • Numerical models support efficient wave-propagation and wave-structure studies, but must represent relevant physics and approximate the mathematical model accurately.
  • Potential-flow assumptions simplify Navier–Stokes modeling but limit applicability to settings such as non-overturning waves and restricted structures.
  • High-order convergent schemes are important for reducing phase lag during long-time wave propagation.
  • Unfitted methods avoid requiring computational meshes to conform to physical interfaces and boundaries, while differing in geometry representation and boundary enforcement.
  • The shifted boundary method keeps intersected elements whole and modifies boundary conditions using boundary mappings and Taylor expansions to approximate true conditions.
  • FNPF wave propagation and wave-structure interaction remain relatively underexplored with unfitted discretizations despite their moving and deforming complex boundaries.
  • The proposed model combines unfitted spectral elements, polynomial-corrected shifted boundaries, polynomial-preserving gradient recovery, arbitrary-order free-surface differences, and hyperviscosity.
  • Verification and validation cover periodic and finite tanks, convergence, computational speed, gradient recovery, long-time nonlinear waves, bathymetry changes, and vertical-wall interactions.

2 Mathematical model

The mathematical model describes a time-dependent two-dimensional fluid domain with a moving free surface, bathymetry, walls, and periodic boundaries. It uses potential flow, a Laplace problem, nonlinear free-surface conditions, pressure-based forces, and relaxation-zone wave initialization or control.

  • The fluid domain is two-dimensional and time-dependent, with free surface, bathymetry, walls, and periodic boundaries enclosing the fluid.
  • The potential-flow approximation assumes incompressible, inviscid, irrotational flow and represents velocity as the gradient of a velocity potential.
  • Governing equations: The velocity potential is obtained from a Laplace problem with prescribed free-surface potential, no-flux bathymetry and wall conditions, and periodic potential boundaries.
  • Free-surface conditions: The vertical free-surface velocity wη = ∂zϕ|z=η links the Laplace solution to the coupled kinematic and dynamic free-surface conditions.
  • Forces: The model computes local pressure from Bernoulli’s equation and integrates pressure over wetted boundaries to obtain instantaneous forces.
  • Wave configurations: Periodic constant-depth cases use stream-function initialization, while finite-tank cases start from zero and solitary waves use spectrally accurate interpolated solutions.
  • Wave generation and absorption: Relaxation zones generate and absorb waves by blending free-surface variables with target values through a high-order shape function and gradual ramping.

3 Numerical discretization: The fluid domain

The fluid domain is discretized on an unfitted Cartesian mesh using a polynomial-corrected shifted-boundary method and high-order spectral elements. Boundary data are transferred between true and surrogate boundaries through mappings and high-order corrections, while the Laplace problem is solved variationally.

  • Unfitted domain construction: The unfitted surrogate domain retains whole background elements instead of cutting boundary-intersecting elements.The surrogate is constructed from a regular quadrilateral background tessellation by identifying interior and boundary-intersecting elements.
  • Shifted-boundary formulation: The shifted-boundary method solves on the surrogate domain and modifies boundary conditions to approximate conditions on the true boundary.A mapping between the true and surrogate boundaries transfers Dirichlet and Neumann data.
  • Boundary mappings: For the free surface, the boundary mapping is purely vertical, whereas walls and bathymetry use mappings along the true normal direction.This mapping freedom accommodates the distinct geometries of free surfaces and solid or bathymetric boundaries.
  • Boundary correction: Naively equating boundary data on the true and surrogate boundaries limits the shifted-boundary treatment to first-order accuracy.The method therefore introduces a more elaborate correction to preserve higher-order accuracy.
  • Boundary correction: The polynomial-corrected formulation transfers boundary data through mapped solution evaluations and surrogate-normal decompositions within a variational penalty framework.The formulation separates normal and tangential contributions and retains penalty stabilization through parameter τ.
  • Spectral-element discretization: The fluid variables use tensor-product spectral elements of polynomial degree p on Cartesian quadrilaterals, producing GLL nodes and a sparse Galerkin linear system.The assembled system has the form Aϕh = b and is solved directly using Gaussian elimination in MATLAB.
  • Resolution choice: At ε/εmax = 90%, the highly nonlinear wave configuration fits three vertical elements, consistent with reported FNPF-SEM refinement levels.The discretization is motivated by resolving vertical dispersion with relatively few elements and moderate-to-high polynomial order.

4 Numerical discretization: The free surface

The free surface is evolved with strong finite differences and explicit Runge–Kutta time integration, while polynomial-preserving gradient recovery computes the vertical velocity and hyperviscosity stabilizes nonlinear simulations.

  • Free-surface discretization: The nonlinear free-surface equations use finite-difference operators on GLL points and are integrated with an explicit fourth-order, four-stage Runge–Kutta scheme.Unless otherwise specified, the constant timestep is Δt = T/200.
  • Gradient recovery: The vertical free-surface velocity wη = ∂zϕ|z=η links the Laplace problem to both free-surface conditions and strongly affects stability.Accurate recovery of this derivative is therefore central to the overall numerical scheme.
  • Gradient recovery: Local recovery is efficient but can lose accuracy and introduce discontinuities across element interfaces, especially when extrapolating to an unfitted boundary.These limitations motivate a polynomial-preserving recovery procedure.
  • Polynomial-preserving recovery: The polynomial-preserving approach fits a one-dimensional vertical least-squares polynomial of order q = p + 1 using the closest q + 1 vertical samples.The coefficients minimize the squared residual between polynomial values and sampled velocity-potential values.
  • Polynomial-preserving recovery: Differentiating the fitted polynomial and evaluating at the free-surface point yields the gradient from coefficient c1, with expected error O(h^q=p+1).The approach incorporates projection operations on unfitted domains and extends to three dimensions with additional free-surface points.
  • Stabilization: Hyperviscosity is added because highly nonlinear waves can develop aliasing-driven instabilities during long-time integration.Its parameters include grid spacing Δx, scale coefficient γ, order m, and a finite-difference free-surface Laplacian Dxx.
  • Stabilization: Unfitted geometrical errors can generate small-scale effects even when numerical and geometric results converge, making stabilization choices consequential.The stabilization parameters are tested for long-time simulations of very nonlinear waves.

5 Results: Verification in a periodic wave flume

Periodic-flume tests examine convergence, long-time accuracy, computational efficiency, gradient recovery, and hyperviscosity stabilization for highly nonlinear waves. The results support high-order approximations, PPR gradient recovery, and stabilization parameters selected from broad long-time tests.

  • Convergence under mesh-refinement: All polynomial orders converge at least as O(h^p), with most approaching O(h^(p+1)) for p ∈{1, 2, 3}.The convergence study covers kh ∈{1, π, 2π} and polynomial orders p ∈{1, . . . , 5}.
  • High-order approximations and gradient recovery: The coarsest mesh is selected to satisfy erel(100T) ≤10%, using PPR and local recovery across polynomial orders p ∈{1, . . . , 6}.The estimate uses one-period convergence behavior and assumes erel(100T) ≈100 erel(1T).
  • High-order approximations and gradient recovery: PPR exhibits roughly O(h^(p+1)) convergence versus O(h^p) for local recovery, although both methods fall slightly below expectations for p ≥4.The deviation is attributed to the fourth-order ERK4 scheme and fixed time-step size.
  • High-order approximations and gradient recovery: A p = 6 PPR run achieved 3.03% error after 100 wave periods in roughly 83 seconds, compared with 9.70% and around 13,319 seconds for p = 2 local recovery, a 1,593 speed-up.These CPU measurements cover only time integration.
  • Long-time simulation with and without stabilization: The tested hyperviscosity settings produced 487 finished runs out of 1,620, while 908 exploded and 225 were too damped.Among the 50 best runs, (m, γ) = (3, 1) occurred 13 times; the study then used (m, γ) = (3, 1).
  • Long-time simulation with and without stabilization: With (m, γ) = (3, 1), diffusive errors were 1–3 orders of magnitude lower than dispersive errors and remained nearly constant over time.Total errors increased with wave nonlinearity for kh ∈{1, π, 2π}.

6 Results: Numerical experiments in a finite wave flume

Finite-flume experiments assess long-time propagation, harmonic generation over a submerged bar, and solitary-wave interactions with walls and curved bathymetry. The unfitted framework shows agreement with reference results across these settings, while extreme wave-wall interaction exposes a breakdown boundary.

  • Experiment scope: The finite-flume studies cover generation and absorption, submerged-bar harmonic transformation, vertical-wall reflection, and propagation over a semi-circular bathymetry bump.Simulations use the unfitted fixed-discretization approach across mild to highly nonlinear wave conditions.
  • Wave generation and absorption: 200 wave periods with kh = 1 and ε/εmax = 70% show minimal profile differences between 100 and 200 periods, indicating effective absorption.The computed profile at x = 7L has reasonable agreement with the stream-function reference, although a slight phase lag is visible.
  • Submerged bar: Submerged-bar propagation steepens the wave on the incline and generates higher-order harmonics, with qualitative agreement with experimental results.Spectral components across the domain support the expected harmonic release and transformed wave field.
  • Vertical wall: Maximum run-up, attachment and detachment heights, and maximum horizontal force show good agreement for solitary waves reflecting from a vertical wall.The inelastic interaction becomes visually apparent for a/h > 0.2.
  • Vertical wall: For a/h = 0.7, the simulation breaks down after predicting only attachment height, the initial force profile, and maximum force.The authors associate the instability primarily with non-uniqueness of the free-surface elevation near the wall under extreme interaction; stabilization might slightly improve results for a/h > 0.6.
  • Curved bathymetry: A mildly nonlinear solitary wave over a semi-circular bump produces total horizontal dynamic-force results that qualitatively agree with digitized numerical FNPF results.The force is obtained by projecting pressure onto the true boundary and integrating along artificial boundary-fitted curves.

7 Conclusion

The conclusion presents an unfitted high-order framework for highly nonlinear water waves and identifies geometric adaptivity as an ongoing requirement for future structure-interaction problems.

  • Conclusion: The framework combines polynomial-corrected shifted-boundary approximation, spectral elements, polynomial-preserving gradient recovery, and wave stabilization.The method targets highly nonlinear waves affected by bathymetric changes.
  • Conclusion: Polynomial-preserving gradient recovery is used to achieve optimal convergence when computing the free-surface vertical velocity.
  • Future work: Future extensions target partially or fully submerged fixed, forced, and freely floating or oscillating structures.The authors emphasize geometric flexibility because these problems contain multiple length scales.
  • Future work: Combining the unfitted methodology with non-conformal element discretization is proposed as a route toward geometric adaptivity on regular Cartesian grids.

Funding

The work contributes to JV’s PhD project on advanced simulation techniques for wave energy converters, funded by DTU Compute.

  • Funding: The PhD project “New Advanced Simulation Techniques for Wave Energy Converts” is fully funded by DTU Compute.

CRediT author contribution statement

All authors contributed to conceptualization, methodology, writing, and visualization, with additional individual roles spanning software, validation, analysis, investigation, and supervision.

  • Contributions: All authors shared conceptualization, methodology, drafting, review and editing, and visualization responsibilities.
  • Contributions: JV additionally handled software, validation, formal analysis, and investigation, while APEK, HBB, and MR provided supervision.

A Comparing effect of stabilization parameters

Table 3 compares stabilization parameters across wave conditions, resolutions, polynomial settings, and nonlinear simulations, ranking the 50 best runs by combined error.

  • Parameter comparison: Table 3 tests stabilization parameters for different waves and resolutions using m, γ, kh, Nx, Nz, and p combinations.The tested sets are m = {1, . . . , 5}, γ = {0, 0.1, 1, 10}, kh = {1, π, 2π}, Nx = {4, 8, 16}, Nz = {1, 2, 4}, and p = {4, 5, 6}.
  • Evaluation condition: The comparison evaluates simulations at ε/εmax = 90% nonlinearity after 50T.This defines the nonlinear and temporal evaluation condition for the analysis in Section 5.3.
  • Ranking criterion: The 50 best runs are sorted using emix = edif + edis.The ranking combines the identified error components edif and edis.

B Free surface wave elevations for the submerged bar test

Figure 13 presents free-surface wave-elevation time series from the submerged bar test at ten wave-gauge locations.

  • Wave-elevation measurements: The figure shows time series of free-surface wave elevation for the submerged bar test.The test is identified as the case discussed in Section 6.2.
  • Gauge locations: Wave elevation is recorded at 10 different wave-gauge locations.The multiple gauge locations support comparison of the simulated wave signal across the test domain.
  • Test configuration: The plotted quantity is organized as a time series for the submerged bar configuration.The figure focuses on wave propagation measurements rather than a static spatial field.

C Unfitted mesh for the semi-circular bathymetry bump

Figure 14 provides a zoomed view of the unfitted mesh around the semi-circular bathymetry bump.

  • Bathymetry geometry: The figure focuses on the unfitted mesh near a semi-circular bathymetry bump.The geometry shown is a curved bathymetric feature represented within the mesh.
  • Local mesh view: Figure 14 is a zoomed-in view of the mesh region surrounding the bump.The zoom isolates the local mesh treatment around the curved feature.
  • Mesh representation: The displayed mesh is explicitly identified as unfitted.This distinguishes the mesh representation used around the semi-circular bathymetry geometry.
Loading 2608.28123v1…