PaperlayerView exact revision on arXiv

Source-linked AI summary

State-of-the-art SPH solver DualSPHysics: from fluid dynamics to multiphysics problems

Jose M. Domínguez, Georgios Fourtakas, Corrado Altomare, Ricardo B. Canelas, Angelo Tafuni, Orlando García-Feal, Ivan Martínez-Estévez, Athanasios Mokos, Renato Vacondio, Alejandro J. C. Crespo, Benedict D. Rogers, Peter K. Stansby, Moncho Gómez-Gesteira

arXiv:2104.00537v1physics.flu-dynmath.NA

TL;DR

DualSPHysics addresses the computational demands of weakly compressible SPH while extending it to engineering-relevant multiphysics and coupled simulations. The paper presents these developments and their validation, reporting increasing device performance with problem size and continued applicability across wave, multiphase, and fluid-driven-body problems. It also identifies remaining limitations associated with truncated kernel support and the weakly compressible formulation.

  • Problem

    Weakly compressible SPH is computationally demanding because of numerous neighbouring-particle interactions and Courant-restricted timesteps.

  • Method

    The paper presents DualSPHysics developments, validation practices, GPU-enabled implementation, multiphase formulations, and coupling with wave, DEM, Chrono, and MoorDyn models.

  • Results

    GPIPS rises from about 0.14 for CPU2 to 29 for the Tesla V100 as particle count becomes large enough to use the hardware efficiently.

  • Takeaways & Limitations

    DualSPHysics supports applications including hydrodynamic loading, natural-hazard risk, floating-body and exposed-structure design, and multiphase sediment-flow modelling.

  • Takeaways & Limitations

    Density-diffusion terms are not consistent near truncated kernel support such as free surfaces or wall boundaries, although one improves wall-pressure behaviour.

Abstract

from arXiv · show

DualSPHysics is a weakly compressible smoothed particle hydrodynamics (SPH) Navier-Stokes solver initially conceived to deal with coastal engineering problems, especially those related to wave impact with coastal structures. Since the first release back in 2011, DualSPHysics has shown to be robust and accurate for simulating extreme wave events along with a continuous improvement in efficiency thanks to the exploitation of hardware such as graphics processing units (GPUs) for scientific computing or the coupling with wave propagating models such as SWASH and OceanWave3D. Numerous additional functionalities have also been included in the DualSPHysics package over the last few years which allow the simulation of fluid-driven objects. The use of the discrete element method (DEM) has allowed the solver to simulate the interaction among different bodies (sliding rocks, for example), which provides a unique tool to analyse debris flows. In addition, the recent coupling with other solvers like Project Chrono or MoorDyn has been a milestone in the development of the solver. Project Chrono allows the simulation of articulated structures with joints, hinges, sliders and springs and MoorDyn allows simulating moored structures. Both functionalities make DualSPHysics one of the meshless model world leaders in the simulation of offshore energy harvesting devices. Lately, the present state of maturity of the solver goes beyond single phase simulations, allowing multi-phase simulations with gas-liquid and a combination of Newtonian and non-Newtonian models expanding further the capabilities and range of applications for the DualSPHysics solver. These advances and functionalities make DualSPHysics a state-of-the-art meshless solver with emphasis on free-surface flow modelling.

1. Introduction

SPH supports violent free-surface and multiphysics simulations but is computationally demanding because of extensive particle interactions and Courant-limited timesteps. DualSPHysics addresses these challenges through open-source GPU acceleration, dual CPU/GPU execution, accessible tooling, and expanded physical capabilities.

  • SPH enables simulations spanning violent coastal and offshore flows, multiphase industrial processes, deformable solids, fluid-structure interaction, and other phenomena.
  • Large numbers of neighbouring particle interactions and Courant-restricted timesteps make weakly compressible SPH computationally demanding.GPU acceleration is particularly suitable because SPH uses vector lists of particles and their interactions.
  • DualSPHysics is an open-source package developed to exploit modern GPUs and enable real engineering simulations on a standard desktop PC.The project reports 52,000 total downloads and up to 10,000 downloads for each sub-release.
  • DualSPHysics runs on CPUs with OpenMP or Nvidia GPUs with CUDA, with CPU features validated before being ported to CUDA.This dual-development strategy is intended to facilitate implementation of new features by developers.
  • Dedicated pre- and post-processing tools support input generation, visualization, and extraction of quantities such as vorticity, surface elevation, and forces.DesignSPHysics, a FreeCAD GUI plugin, makes generating new test cases more straightforward and accessible.
  • Version 5.0 added DEM and Project Chrono coupling, MoorDyn coupling, multiphase solvers, wave-generation features, wave-model coupling, and inlet/outlet conditions.The release also includes specialized second-order wave generation, active absorption, and relaxation zones.
  • The paper presents DualSPHysics developments and emphasizes validation before software releases, covering formulations, features, couplings, implementation, applications, and challenges.

2. SPH formulation

The SPH formulation represents smooth fields through compactly supported kernel convolutions and their discrete particle approximations. DualSPHysics implements multiple kernel choices, including cubic B-spline and Wendland C2 functions, with particle sets covering fluid, boundary, and floating-object particles.

  • 2.1. Interpolants and kernel functions: SPH represents a smooth function over a domain as a convolution with a kernel function.
  • 2.1. Interpolants and kernel functions: A kernel is a smooth, non-negative, normalized function of q = |r|/h with compact support, vanishing when q reaches the support radius k.
  • 2.1. Interpolants and kernel functions: The discrete SPH approximation replaces the continuous convolution with a sum over neighbouring particles and associated particle volumes.
  • 2.1. Interpolants and kernel functions: The particle set P contains fluid particles F, boundary particles B, and floating objects K, with P = F ∪ B and K ⊂ B.
  • 2.1. Interpolants and kernel functions: SPH gradients are obtained through analogous continuous and discrete kernel-gradient expressions based on convolution identities.
  • 2.1. Interpolants and kernel functions: The kernel approaches the Dirac delta function as h tends to zero, and several kernel functions are implemented in DualSPHysics.
  • 2.1. Interpolants and kernel functions: DualSPHysics implements the third-order B-spline cubic kernel, whose normalization coefficient α_D is 10/7πh^2 in 2-D and 1/πh^3 in 3-D.
  • 2.1. Interpolants and kernel functions: The Wendland C2 kernel has a positive Fourier transform and better pairing-instability characteristics, with α_D equal to 7/4πh^2 in 2-D and 21/16πh^3 in 3-D.

2.2. Governing equations

DualSPHysics formulates compressible-fluid dynamics with Lagrangian Navier-Stokes equations. The continuity and momentum equations describe density evolution and velocity evolution under pressure, dissipation, and external forces.

  • The governing equations are the Navier-Stokes equations for a compressible fluid written in Lagrangian form.
  • The momentum equation evolves velocity using pressure gradients, dissipation terms, and accelerations from external forces such as gravity.

2.3. SPH discretisation of the governing equations

DualSPHysics discretizes conservation laws with particle-based SPH operators for density, momentum, diffusion, viscosity, and sub-particle-scale turbulence. These formulations include stabilizing terms but have boundary-support limitations.

  • SPH conserves particle mass and computes density changes through the discrete continuity equation, including a numerical density-diffusion term.
  • The discrete Euler equation uses a symmetric SPH operator that guarantees momentum conservation.
  • Two density-diffusion formulations filter high-frequency numerical noise, smoothing density and pressure to improve weakly compressible SPH stability.
  • Both density-diffusion terms are inconsistent near truncated kernel support, although the dynamic-pressure formulation improves pressure behavior near walls.
  • Artificial viscosity reduces oscillations and stabilizes the scheme, while laminar viscous dissipation is approximated using the kinematic viscosity.
  • Sub-particle scale model: The implemented SPS turbulence model uses an eddy-viscosity closure with Smagorinsky constant C_S=0.12 and C_l=0.0066.

2.4. Equation of state and compressibility

DualSPHysics couples density and pressure through an equation of state for weakly compressible flow. The numerical sound speed is selected to limit density variation while permitting larger explicit time steps.

  • Density and pressure are coupled through an equation of state based on the numerical speed of sound, enabling weak compressibility.
  • c_s=10‖v‖_max targets 1% density variation and Ma≈0.1, allowing larger time steps in explicit integration.

2.5. Time integrators and time step

DualSPHysics provides Verlet-based explicit integrators for advancing particle variables. The CFL condition bounds the variable time step to maintain the physical domain of dependence.

  • DualSPHysics implements two explicit time-integration schemes for advancing density, position, velocity, and acceleration.
  • Verlet scheme: The Verlet scheme is low-cost, second-order accurate in space, and avoids multiple calculation steps within an iteration interval.
  • Verlet scheme: Because staggered integration decouples density and velocity, an intermediate step is recommended every N_S≈40 steps to prevent divergence.
  • Symplectic position Verlet scheme: The symplectic position Verlet scheme is second-order accurate in time, time reversible, and symmetric without diffusive terms.
  • Time step: The CFL condition bounds explicit integration, with the variable time step selected as the minimum of the relevant limits and C_CFL usually between 0.1 and 0.2.

2.6. Boundary conditions

DualSPHysics represents solid boundaries with particles that can move according to prescribed or interaction-driven forces. Its mDBC formulation introduces ghost nodes to address limitations of the dynamic boundary condition and impose slip choices.

  • Solid boundaries use boundary particles distinct from fluid particles, and their positions can be moved according to user-defined motion.
  • Dynamic boundary condition: The dynamic boundary condition evaluates continuity on wall particles and converts density increases into repulsive fluid-boundary pressure forces.
  • Modified dynamic boundary condition: From version 5.0, mDBC projects a ghost node across the boundary, computes fluid properties there with corrected SPH, and mirrors them to boundary particles.
  • Modified dynamic boundary condition: mDBC supports no-slip boundaries by reversing ghost-node velocity and free-slip boundaries by assigning the exact tangential velocity.

2.7. Periodic boundary conditions

DualSPHysics implements periodic boundary conditions by cyclically connecting particles near opposite open lateral boundaries. It also uses Fickian particle shifting to reduce discretisation errors from anisotropic particle distributions.

  • 2.7. Periodic boundary conditions: Periodic boundary conditions represent an infinitely long domain by connecting particles within 2h of opposite lateral boundaries in 2-D and 3-D.The smoothing kernels use k=2, so particles near one open lateral boundary interact with particles across the opposite boundary.
  • 2.7. Periodic boundary conditions: Particle shifting maintains a near-isotropic distribution to reduce discretisation errors from zeroth- and higher-order kernel moments.The Fickian-based algorithm is used especially when particles become anisotropically distributed during violent flows.
  • 2.7. Periodic boundary conditions: The shifting distance is computed from the kernel gradient and a diffusion coefficient proportional to particle velocity magnitude.The method restricts particle movement to prevent excessive displacement and loss of information between domain cells.

3. Functionalities

DualSPHysics extends its SPH formulation with wave-generation and absorption tools, fluid-driven rigid-body dynamics, versatile open boundaries, and multiphase non-Newtonian capabilities. These features broaden the solver from free-surface flow modelling toward coupled engineering applications.

  • Wave generation and absorption: Moving boundary particles generate regular, random long-crested, and solitary waves using piston- and flap-type wavemakers.Active and passive absorption techniques are included, with active control correcting wavemaker motion to reduce reflections.
  • Wave generation and absorption: Second-order wavemaker formulations and spectrum-based methods support regular, irregular, and multiple solitary-wave generation.Implemented irregular-wave spectra include JONSWAP and Pierson-Moskowitz; second-order methods target parasitic and spurious long waves.
  • Fluid-driven objects: Fluid-driven objects move through forces from surrounding fluid particles followed by rigid-body integration of translational and angular motion.The resulting particle velocity combines object translation with angular velocity crossed with the offset from the centre of mass.
  • Open boundaries: Open-boundary regions can operate simultaneously as inlets and outlets, supporting backflow, flow reversion, and unsteady velocity or pressure profiles.The implementation also supports variable free-surface elevation, which is important for flows with changing water levels.
  • Multiphase modelling: The multiphase formulation addresses liquid-gas flows and supports combinations of viscous operators, gradients, Newtonian models, and non-Newtonian models.A selectively activated free-surface correction adapts particle shifting to prevent unnatural voids and maintain stability for high-density-ratio violent flows.

4. Coupling with other models

DualSPHysics couples with wave-propagation, DEM, Project Chrono, and MoorDyn models to combine local SPH resolution with external descriptions of waves, contacts, articulated bodies, and moorings. The reported validation results support coupling accuracy at approximately the smoothing-length scale in 2-D tests.

  • Coupling rationale: Coupling SPH with other numerical methods is intended to extend both models beyond their individual application fields for wider-ranging simulations.The motivation includes reducing the computational burden of complex, long-duration, multiscale coastal simulations.
  • Coupling with wave propagation models: SWASH and OceanWave3D coupling schemes use piston, relaxation-zone, or open-boundary interfaces to transfer wave velocities and surface elevations.The OceanWave3D approach is two-way: DualSPHysics sends surface elevation back after imposing interface data and applying wave-reflection correction.
  • Validation: Validation of the coupling methodologies used 2-D tests and showed model accuracy generally on the order of the smoothing length.Open boundaries enable shallow-water coupling where wave nonlinearity and mass transport dominate.
  • Coupling with DEM: DCDEM represents simultaneous solid-fluid and solid-solid interactions for granular-type engineering flows.Applications include hydrodynamic loading, natural-hazard risk, and floating- or exposed-structure design, subject to the approach’s limitations.
  • Coupling with Project Chrono: Project Chrono solves constrained multibody dynamics involving contacts, joints, hinges, sliders, and springs, while DualSPHysics updates floating particles from Chrono’s results.DualSPHysics passes accelerations and the timestep to Chrono, which returns constrained positions and velocities.
  • Coupling with MoorDyn: MoorDyn receives fairlead kinematics, computes mooring-line tensions, and transfers those tensions back as external forces in DualSPHysics.This coupling was applied to the motion of a moored floating structure under waves.

5. Code and implementation

DualSPHysics provides CPU and GPU implementations organized around neighbour-list, particle-interaction, and system-update tasks. Its performance evaluation demonstrates substantial GPU scaling for large particle counts, while the open-source codebase supports compilation and extension across platforms.

  • 5.1. Implementation: DualSPHysics runs on CPU and GPU architectures, with the same formulation optimized for both and CUDA used for Nvidia GPUs.CPU execution uses C++, while GPU execution uses CUDA kernels and can exploit large GPU computational capacity.
  • 5.1. Implementation: Each simulation step repeats neighbour-list construction, particle interaction, and system update tasks until the target physical time is reached.The system update advances particle position, density, and velocity using current properties, interaction results, and timestep duration.
  • 5.1. Implementation: Particle interaction exceeds 90% of total execution time because it solves momentum and continuity equations through neighbour interactions.Cell-linked lists limit each particle’s interactions to neighbouring particles in adjacent cells.
  • 5.2. Source files and compilation: The code contains more than 240 source files and over 95,000 lines, organized into core, functionalities, I/O, and general-purpose blocks.Included build systems support Linux, Windows, and CMake compilation, including compilation without CUDA when no GPU is available.
  • 5.3. Performance: 80,353,962 particles were simulated in under 48 hours on an NVIDIA Tesla V100, compared with 2.15 million on the fastest tested CPU.The evaluation used the second SPHERIC validation dam-break test over the first two physical seconds.
  • 5.3. Performance: 29 GPIPS was achieved by the Tesla V100, and its speedup against CPU2 exceeded 200x from 4 million particles onward.The Tesla V100 reached saturation at 12 million particles, while the reported comparison against the RTX 2080 Ti was about 1.3x.

[MHz]

DualSPHysics combines hardware-accelerated simulation with integrated tools for case generation, execution, visualization, and analysis. Its workflow supports accessible model construction through DesignSPHysics.

  • Performance: Using cells of size h improves neighbour search, enabling CPU and GPU speedups of 1.14x and 1.3x, respectively.This optimisation increases memory consumption and reduces the maximum number of particles simulated on a GPU.
  • Performance: Using cells of size 2h supports 6.05 million particles per Gigabyte of device memory, compared with 5.8 million particles per Gigabyte for h.GPU memory size remains a major limitation for simulations requiring large particle counts.
  • Performance: A multi-GPU MPI approach removes the GPU-memory limitation and enables simulations with more than 1,000 million particles at relatively short runtimes.The approach combines GPUs in the same or different computation nodes.
  • Workflow: DualSPHysics provides a standalone workflow spanning initial-condition generation, SPH solving, and physical-magnitude analysis and visualization.The package includes pre-processing and post-processing tools alongside the solver.
  • Workflow: GenCase deploys particles on a 3-D Cartesian mesh, while XML inputs specify geometries, constants, and solver options.Changing the initial particle distance allows the same geometry to be discretized at different resolutions.
  • Workflow: DesignSPHysics integrates pre-processing, solving, and post-processing in a FreeCAD GUI intended to help new users create and run cases from scratch.It is open-source and uses Python and QT within FreeCAD.

6. Applications

DualSPHysics has been validated across coastal, hydraulic, fluid-structure, wave-energy, and multiphase applications. Its extensions couple SPH with propagation, DEM, mooring, rigid-body, and non-Newtonian models.

  • Validation and application scope: DualSPHysics features have been validated against analytical solutions, numerical benchmarks, and experimental results, with Table 4 collecting recent validated applications.The literature includes coastal structures, multiphase flows, fluid-driven objects, wave-energy converters, and hydraulic problems.
  • Coastal engineering: Coastal applications include validated second-order wave generation and absorption, porous-breakwater run-up, real sea-dike analysis, and coupled wave propagation.OceanWave3D coupling reduces the cost of propagating incoming waves far from coastal structures.
  • Fluid-driven objects: DEM coupling enables simulations of solid-solid interactions in free-surface flows, including dam-break impacts and debris-flow problems.The coupled framework solves fluid flow and rigid-body dynamics within the same meshless framework and was experimentally validated.
  • Fluid-driven objects: MoorDyn coupling reproduced moored-box motions and line tensions in good agreement with scale-model tests across different regular waves.Further validation is still needed for several sea states, including irregular waves and wave trains near the structure’s natural frequency.
  • Wave-energy converters: Project Chrono coupling supports articulated wave-energy devices, and an OWSC simulation accurately predicted hydrodynamic response under regular and irregular waves.The model included mechanical constraints, hydraulic PTO nonlinearities, and frictional contacts between the flap and bearings.
  • Hydraulic engineering: A 3-D vertical-slot-fishway simulation using buffer-region open boundaries agreed well with experiments for discharges, velocity profiles, and water elevations.The application represents free-surface turbulent subcritical flow.
  • Multiphase modelling: Multiphase capability addresses liquid-gas and liquid-sediment interactions, while a shifting algorithm was developed to handle unphysical voids formed during accelerated high-resolution simulations.Applications include dam breaks, fuel-tank sloshing, and erosion-related non-Newtonian sediment flow.

7. Outlook

The paper presents DualSPHysics as a mature open-source SPH code whose developments target boundary conditions, solver coupling, and industrial applicability. Its outlook is constrained by the computational demands of the resolution needed for engineering physics.

  • Outlook: DualSPHysics v5.0 consolidates validated developments in SPH formulation, code structure, coupling, and applications summarized across the package.The developments address SPHERIC Grand Challenges GC2-Boundary conditions, GC4-Coupling, and GC5-Applicability to industry.
  • Outlook: The paper identifies industrial uptake of SPH as requiring special efforts, particularly because industrial applicability has received less attention than it deserves.The conclusion frames broader use in industry as an ongoing development priority.
  • Outlook: Engineering simulations may require resolutions beyond the computational capabilities of GPU-accelerated or massively parallel CFD software packages.SPH’s Lagrangian formulation can make the resolution and computational demands of implementable physical processes especially challenging.
  • Outlook: DualSPHysics is highly optimized for Nvidia GPU hardware through CUDA kernels, templates, and other optimizations.The implementation reflects continuing changes in GPU and CPU computing architectures.
Loading 2104.00537v1…