A high-order Fourier Continuation (FC)-based spectral incompressible Smoothed Particle Hydrodynamics (ISPH) scheme for general boundary conditions in wall-bounded domains

This paper introduces a high-order Fourier Continuation (FC)-based spectral incompressible Smoothed Particle Hydrodynamics (ISPH) scheme that extends the method to wall-bounded domains with general boundary conditions, enabling high-order convergence and accurate simulation of complex vortex dynamics through frequency-space discretization on a periodic extension of the domain.

The Gist

Imagine you are trying to simulate how water flows through a pipe or how a swirl of smoke moves near a wall. To do this on a computer, scientists use a method called Smoothed Particle Hydrodynamics (SPH). Think of SPH as a digital crowd of tiny, invisible marbles. Instead of using a fixed grid (like graph paper), the computer tracks these marbles as they move, bounce, and swirl around.

For a long time, there was a problem with a specific, super-accurate version of this method called "Spectral SPH." It was like having a super-fast sports car that could only drive on a perfectly circular racetrack. If you tried to drive it on a straight road with walls (like a pipe), the math would break down, creating "ghosts" or glitches in the simulation. This is because the math behind this method loves periodicity—it assumes the world loops around forever, like a Pac-Man screen where if you go off the right edge, you appear on the left.

But real life isn't like Pac-Man. Real pipes have walls where water stops or slides, and smoke doesn't loop around the room.

The Solution: The "Magic Extension" (Fourier Continuation)

The authors of this paper, Meixuan Lin and colleagues from the University of Manchester, invented a clever trick called Fourier Continuation (FC) to fix this.

Here is the analogy:
Imagine you are trying to sing a song that loops perfectly, but you have a verse that ends abruptly on a high note. If you try to loop it, it sounds like a jarring screech.

• The Old Way: You just cut the song and loop it. It sounds terrible (mathematically, this is called the "Gibbs phenomenon").
• The New Way (FC): Before you loop the song, you add a short, smooth "bridge" at the end. You write a few extra notes that gently bring the high note down to match the starting note, creating a seamless loop.

In the computer simulation, the researchers do this mathematically:

1. Fitting: They look at the data right next to the wall (the "end of the song").
2. Extrapolating: They use a high-order polynomial (a fancy mathematical curve) to predict what the data would look like if it continued past the wall.
3. Blending: They smoothly mix this prediction with the data on the other side of the wall to create a seamless, smooth loop.

By doing this, they trick the computer into thinking the wall is just part of a giant, smooth, looping world. This allows the "super-fast sports car" (the spectral method) to drive on the straight road (the wall-bounded domain) without crashing.

What They Tested

To prove their "magic extension" works, they ran several tests:

• The Gaussian Vortex: They simulated a perfect swirl of wind moving across the screen. Without their trick, the swirl would get distorted when it hit the edge. With the trick, it flowed smoothly off the screen, just like a real wind gust.
• Poiseuille Flow: This is water flowing through a pipe, pushed by a constant force. The math for this is a simple curve. Their method predicted this curve with incredible precision, better than standard methods.
• Couette Flow: Imagine two parallel plates, one stationary and one moving, with fluid in between. The fluid has to match the speed of the moving plate and stop at the stationary one. This is a tricky "asymmetric" problem. Their method handled it naturally, without needing complex workarounds.
• The Vortex Rebound: This is the "boss level" test. They simulated two spinning whirlpools crashing into a wall. When they hit, they create tiny, complex secondary swirls and bounce back. This is very hard to simulate accurately. Their method matched the results of other top-tier, highly accurate scientific software, proving it could capture these tiny, complex details.

The Result

The paper concludes that by adding this "magic extension" (Fourier Continuation), they have successfully upgraded the Spectral SPH method.

• Speed: It remains very fast (using a mathematical shortcut called FFT).
• Accuracy: It is "high-order," meaning it gets much more precise as you add more particles, capturing fine details like tiny vortices.
• Versatility: It can now handle walls, inflows, and outflows, not just looping worlds.

The Catch (Limitations)

The authors are honest about the current limits. Right now, this "magic extension" works best on simple, smooth, rectangular shapes (like a straight pipe or a box). It doesn't yet work well on complex, jagged shapes like a car engine or a human heart. They plan to fix this in future work to make it a truly universal tool for any shape.

In short, they found a way to make a super-accurate, fast fluid simulation method work in the real world, where walls exist and things don't loop forever.

Technical Summary

Technical Summary: A High-Order Fourier Continuation (FC)-Based Spectral ISPH Scheme

Problem Statement
Smoothed Particle Hydrodynamics (SPH) is a robust meshfree method well-suited for problems involving large deformations and free surfaces. However, achieving high-order accuracy and computational efficiency remains a significant challenge within the SPH community. While previous work by the authors introduced a spectral SPH scheme that reformulates operators in the frequency domain using the Fast Fourier Transform (FFT) to achieve $O(N \log N)$ complexity, this approach was inherently restricted to periodic domains. The periodicity constraint of the FFT prevents the direct application of the spectral SPH scheme to wall-bounded flows with non-periodic boundary conditions, such as Dirichlet (no-slip) or Neumann conditions. Existing high-order SPH methods typically operate in physical space and face difficulties in matching the spectral accuracy and efficiency of frequency-domain approaches.

Methodology
To overcome the periodicity limitation, this paper introduces a high-order Fourier Continuation (FC) algorithm into the spectral incompressible SPH (ISPH) framework. The core methodology involves the following components:

1. Fourier Continuation (FC) Extension:
   The non-periodic computational domain is extended to create a smooth, periodic function suitable for FFT-based processing. This is achieved through a three-step polynomial-based process:

   • Boundary Polynomial Fitting: High-order polynomials are fitted to grid points near the boundaries using a least-squares approach.
   • Extrapolation: These polynomials are extrapolated into an extension region of length $d$ (defined as a fraction of the domain size).
   • Smooth Blending: A smooth blending function (a 5th-order polynomial) is used to seamlessly link the extrapolated values from opposite boundaries, ensuring the extended function is $C^p$ smooth and periodic.
   • Neumann Condition Handling: For pressure fields satisfying homogeneous Neumann boundary conditions ($\partial P/\partial n = 0$), the polynomial fitting is constrained to enforce a zero derivative at the boundary.

2. Spectral SPH Discretisation:
   The SPH operators are reformulated in the frequency domain. Unlike traditional physical-space SPH which uses linear convolution, this scheme employs circular convolution on the FC-extended domain. This avoids the Gibbs phenomenon associated with zero-padding non-periodic functions. The derivatives are computed via FFT, reducing computational complexity to $O(N \log N)$.

   • Kernels: The study utilizes high-order Gaussian kernels (specifically $G_4$, $G_6$, and $G_8$) to maintain high-order convergence. The resolving power of these kernels is analyzed, showing that higher-order kernels ($G_8$, $G_{10}$) better approximate the Dirac delta function in the frequency domain, extending the range of resolved wavenumbers.

3. Time Integration and Pressure Solver:
   The scheme employs a projection-based time integration method (using a low-storage 3rd-order Runge-Kutta scheme) to decouple velocity and pressure. The incompressibility constraint is enforced by solving a Pressure Poisson Equation (PPE). For wall-bounded flows, the PPE is solved using a combination of Discrete Fourier Transforms (DFT) in the periodic streamwise direction and Discrete Cosine Transforms (DCT-I) in the wall-normal direction, efficiently handling the Neumann boundary conditions on the physical domain.

Key Contributions

• Extension of Spectral ISPH to Wall-Bounded Flows: The primary contribution is the successful integration of the FC algorithm into the spectral SPH framework, enabling the simulation of incompressible flows with arbitrary wall boundary conditions (Dirichlet and Neumann) while retaining high-order spectral accuracy.
• High-Order Convergence: The paper demonstrates that the FC-based spectral SPH scheme achieves convergence rates exceeding 4th order (limited by the kernel order and polynomial fitting degree) for non-periodic domains.
• Efficient Implementation: The method utilizes circular convolution on the extended domain, avoiding the discontinuities that cause Gibbs oscillations in standard spectral methods. The computational cost is shown to be acceptable, with an extension length of $d=0.25N$ adding only a 25% overhead to the domain size while significantly improving accuracy.
• General Boundary Condition Handling: The framework naturally handles asymmetric boundary conditions (e.g., moving walls in Couette flow) without requiring complex splitting methods or homogeneous decompositions, as the FC algorithm smooths any non-periodic function into a periodic one.

Results and Validation
The proposed scheme was validated against several classical CFD benchmarks:

• Convergence Tests: Numerical tests on gradient and Laplacian operators confirmed convergence rates consistent with the polynomial fitting order (5th order) and kernel consistency, with $L_2$ errors decreasing as the extension length increased.
• Gaussian Vortex Convection: The scheme accurately simulated the advection of a Gaussian vortex through a non-periodic domain, demonstrating the ability to handle inflow-outflow conditions without spurious reflections or loss of spectral accuracy.
• Poiseuille and Couette Flows: The method accurately reproduced analytical solutions for both pressure-driven (Poiseuille) and shear-driven (Couette) flows. Notably, the Poiseuille case exhibited "super-convergence" (exceeding 4th order) due to the quadratic nature of the analytical solution.
• Vortex Dipole Rebound: The scheme was tested on the complex vortex-wall interaction problem at Reynolds numbers $Re=100$ and $Re=625$.
  • At $Re=100$, results showed excellent agreement with a high-order Spectral Element Method (SEM) reference (Nektar++), accurately capturing secondary vortex generation and total kinetic energy decay.
  • At $Re=625$, the scheme successfully resolved thin shear layers and secondary vortex detachment. While a resolution of $200 \times 200$ showed some diffusion, a $600 \times 600$ resolution closely matched reference data from Chebyshev pseudospectral methods, accurately capturing peak vorticity magnitudes and locations.

Significance and Claims
The paper claims that the incorporation of Fourier Continuation techniques represents a significant step toward a "fully high-order spectral Lagrangian SPH solver." The work demonstrates that spectral SPH can be extended beyond periodic domains to handle complex wall-bounded flows with high-order accuracy and computational efficiency. The authors emphasize that the method accurately captures complex vortex dynamics and wall interactions, which are critical for high-fidelity fluid simulations.

However, the authors modestly acknowledge current limitations: the scheme is currently restricted to regular, smooth computational domains and relies on Cartesian particle distributions. Consequently, its application to complex geometries is not yet possible. The authors state that future work will address these limitations by adopting a two-dimensional FC framework for general smooth domains, aiming to enable the solver for complex real-world applications.