Smoothed Particle Hydrodynamics Operators

• Smoothed Particle Hydrodynamics (SPH) operators are nonlocal, kernel-based formulations that replace traditional differential operators, enabling meshless simulation of fluid dynamics and continuum mechanics.
• They utilize gradient, divergence, and Laplacian approximations with correction schemes to maintain consistency and robustness even with irregular particle distributions.
• Extensions to SPH operators include nonlocal, fractional, quantum, and anisotropic frameworks, offering versatile solutions for multiphysics, variable-resolution, and high-performance computational challenges.

Smoothed Particle Hydrodynamics (SPH) Operators

Smoothed Particle Hydrodynamics (SPH) is a meshless Lagrangian method for the numerical solution of partial differential equations, primarily those arising in fluid dynamics and continuum mechanics. The central concept in SPH is the replacement of classical differential operators—gradient, divergence, Laplacian—by nonlocal, kernel-weighted integral formulations that can be straightforwardly discretized as sums over particle neighborhoods. The resulting operators form the backbone of SPH simulation methodologies and have been extensively analyzed and generalized for accuracy, stability, and consistency.

1. Nonlocal and Kernel-Based Operator Formulation

The SPH methodology is characterized at the continuum level by the use of nonlocal integral operators, which serve as integral relaxations of local derivatives. For a function $A(\mathbf{x})$ on a domain $\Omega$, the smoothed (kernel-integrated) representation is

$$A(\mathbf{x}) = \int_{\Omega} A(\mathbf{x}')\,W(\mathbf{x}-\mathbf{x}',\,h)\,d\mathbf{x}' + \mathcal O(h^2)$$

where $W(\mathbf{r},h)$ is a nonnegative, normalized smoothing kernel with characteristic width $h$ (Cossins, 2010). The corresponding discrete particle approximation is

$$A_i = \sum_j m_j\,\frac{A_j}{\rho_j}\;W_{ij}\,, \qquad W_{ij} \equiv W(\mathbf{x}_i - \mathbf{x}_j, h)$$

with $m_j$ and $\rho_j$ the mass and density of particle $j$.

Du and Tian introduced rigorous nonlocal versions of the major vector operators for the incompressible Stokes system: a nonlocal diffusion operator $L_\delta[u]$, gradient operator $\nabla_\delta$, and divergence operator $D_\delta$, all parameterized by a nonlocality scale $\delta$. These operators are defined as

$$L_\delta[u](x) = \int_{\mathbb{R}^d} \omega_\delta(|y-x|)[u(y) - u(x)]\,dy$$

$$\nabla_\delta p(x) = \int_{\mathbb{R}^d} (y-x)\,\widetilde\omega_\delta(|y-x|)[p(y) - p(x)]\,dy$$

$$D_\delta u(x) = \int_{\mathbb{R}^d} (y-x)\cdot\widetilde\omega_\delta(|y-x|)[u(y) + u(x)]\,dy$$

with appropriate compactly-supported, radially symmetric kernels $\omega_\delta$ and $\widetilde\omega_\delta$ (Du et al., 2018). These formulations underpin much of the theoretical understanding of SPH operator consistency and nonlocality.

2. Discrete Differential Operators: Gradient, Divergence, and Laplacian

The discrete versions of the SPH differential operators are derived by differentiating the kernel sum or applying vector calculus identities. The standard unsymmetrized and pairwise-symmetric forms are given by (Cossins, 2010, Sutti, 2022):

• Gradient of scalar $A$ (direct, unsymmetric):

$$\nabla A_i = \sum_j m_j\,\frac{A_j}{\rho_j}\,\nabla_i W_{ij}$$

• Gradient (pairwise-symmetric, conservative):

$$\nabla A_i = \rho_i \sum_j m_j \left[\frac{A_j}{\rho_j^2} + \frac{A_i}{\rho_i^2}\right]\nabla_i W_{ij}$$

• Divergence of vector $\mathbf{v}$ (pairwise-symmetric):

$$(\nabla\cdot\mathbf{v})_i = \frac{1}{\rho_i} \sum_j m_j (\mathbf{v}_j - \mathbf{v}_i)\cdot\nabla_i W_{ij}$$

For the Laplacian, two prevalent discretizations dominate:

• Direct second-derivative kernel:

$$\nabla^2 A_i = \sum_j m_j\,\frac{A_j}{\rho_j}\,\nabla^2 W_{ij}$$

• Brookshaw form (first-derivative only):

$$\nabla^2 A_i \approx 2\sum_j m_j\,\frac{A_j - A_i}{\rho_j} \frac{\mathbf{r}_{ij}\cdot\nabla_iW_{ij}}{|\mathbf{r}_{ij}|^2}$$

The Brookshaw form is widely favored for its robustness in the presence of particle noise and is especially prevalent in SPH formulations of viscous and diffusion terms (Cossins, 2010, Lukyanov et al., 2017).

Advanced elliptic SPH operators, notably those required for variable-coefficient problems, have been developed following the two-point-flux and correction frameworks. Schwaiger's correction and the new optimal two-point-flux scheme enforce first-order accuracy throughout the domain and allow for strict monotonicity even under coefficient discontinuities, and are now a standard for stable SPH discretization of elliptic operators (Lukyanov et al., 2017).

3. Consistency, Correction Schemes, and High-Order SPH Operators

Standard SPH operators are only strictly consistent (i.e., reproduce constant, linear, or quadratic fields exactly) for uniform, symmetric particle distributions and constant smoothing length. Several correction schemes have been devised to restore first- or higher-order consistency:

• Gradient correction (e.g., CSPM, Bonet-Lok, moving least squares): Correction matrices are constructed to guarantee reproduction of linear fields by solving local matrix equations, with the corrected kernel gradient

$$\nabla^*_i W_{ij} = \mathbf{L}_i^{-1}\nabla_i W_{ij}$$

where $$\mathbf{L}_i = -\sum_j m_j/\rho_j\,\mathbf{r}_{ij}\otimes\nabla_iW_{ij}$$ (Ghorbani et al., 7 May 2025).

• Consistent second-order gradient and Laplacian: Local correction matrices $G_i$ and $L_i$ ensure exact recovery of gradients on linear fields and Laplacians on quadratic fields, maintaining $O(h^2)$ accuracy even on multiresolution, irregular, or nonuniform particle sets (Hu et al., 2017).
• Multi-dimensional Hessian and anisotropic Laplacian operators: In the anisotropic SPH (ASPH) method, full second-derivative tensors (Hessian matrices) are constructed by suitable kernel and coordinate transformations, supporting accurate simulation of strongly anisotropic diffusion and transport processes (Tang et al., 2024).

A key principle is the satisfaction of normalization/moment conditions on the kernel, ensuring physical consistency and supporting rigorous asymptotics. These include zeroth-, first-, and second-moment constraints. Fulfilling these criteria allows the discrete nonlocal SPH operators to recover their local differential counterparts as the smoothing length $\delta$ or $h$ tends to zero (Du et al., 2018).

4. Operator Extensions: Nonlocal, Fractional, Quantum, and Generalized Frameworks

Nonlocal and Fractional-Order Operators

Nonlocal operator frameworks, as formalized by Du and Tian, provide the continuum foundation for SPH operator consistency, well-posedness, and error estimates (Du et al., 2018). Operators for fractional-order calculus—Riemann–Liouville and Caputo derivatives/integrals—have been formulated in SPH as weighted sums with modified kernel weights carrying explicit fractional and singularity factors. These enable SPH to solve both constant- and variable-order fractional PDEs with high fidelity (Ghorbani et al., 7 May 2025).

Quantum SPH Operators

Quantum algorithms for SPH have demonstrated that all key SPH operator sums—values, gradients, Laplacians—can be encoded as inner products between amplitude-encoded quantum states. Derivative operators correspond to preparing quantum states from kernel derivatives, with exponentially fast convergence in the number of qubits demonstrated for one-dimensional kernels and both Gaussian and Wendland types (Au-Yeung et al., 2020).

Generalized Coordinates and Anisotropy

Total Lagrangian and generalized-coordinate SPH frameworks apply chain-rule operator transformations (via Jacobians) to support nonuniform, coordinate-mapped, or overset domains. Kernels and gradients are defined in reference generalized space $\xi$ and mapped to physical space $x$ via the Jacobian $J$, supporting arbitrary meshless adaptation and high-order accuracy in complex domains (Deng et al., 2021, Tang et al., 2024).

5. Conservation, Stability, and Specializations

Pairwise-symmetric (conservative) forms of the gradient and divergence operators guarantee exact conservation of linear momentum (by antisymmetry), angular momentum (forces along the centers), and, when derived from a particle Lagrangian, total energy (via Noether's theorem) (Cossins, 2010). The rpSPH variant, based on "relative" pressure discretization, enforces Newton's first law by only applying pressure forces when a true gradient is present, eliminating spurious surface tension and particle noise but sacrificing exact momentum conservation—a tradeoff that enables robust handling of instabilities and high Reynold's number flows under certain resolution conditions (Abel, 2010).

Stability and monotonicity in diffusive or elliptic operators are enhanced by two-point-flux formulations that guarantee M-matrix properties, strict positivity of off-diagonal operator weights, and resilience against heterogeneous coefficients and boundary deficiencies (Lukyanov et al., 2017).

Boundary conditions in SPH operator frameworks require particular attention, often using boundary layers, virtual particles, or nonlocal volume constraints to ensure well-posedness and maintain consistency in the presence of truncated kernel support (Du et al., 2018, Hu et al., 2017).

6. Implementation Considerations and Computational Aspects

SPH operator implementations involve the careful choice of kernel, neighbor-search algorithms, data structures for particle properties (mass, position, density), handling of variable smoothing lengths, and time integration schemes (e.g., leapfrog, Runge-Kutta-Fehlberg) (Cossins, 2010, Sutti, 2022). Variable-resolution and particle regularization—via splitting/merging, renormalization, or particle shifting—are essential for maintaining accuracy and stability in large-deformation and multiscale applications (Hu et al., 2017, Deng et al., 2021).

Practical SPH codes, such as the SPHM MATLAB package, employ explicit forms for each operator (density, gradient, divergence, viscous-Laplacian), pairwise-symmetric summations for conservation, and automated neighbor list updates. For elliptic or diffusive problems, linear solvers (e.g., Krylov methods with preconditioning) exploit the structure of the global SPH operator matrices (Sutti, 2022, Lukyanov et al., 2017).

7. Applications and Operator Generalizations

SPH operator frameworks have been extended to:

• Fully anisotropic diffusion, transport, and convection-diffusion systems using full Hessian and tensorial Laplacian discretizations, including in membrane, cardiac, thin-structure, and multiphysics couplings (Tang et al., 2024).
• Variable resolution in multi-domain, fluid-structure, and interface-rich problems via locally-adapted operators and coupling schemes (Hu et al., 2017).
• Fractional-order continuum models for anomalous transport and viscoelastic phenomena (Ghorbani et al., 7 May 2025).
• Quantum and hybrid quantum-classical methods for exponential performance improvements in operator evaluation (Au-Yeung et al., 2020).
• Eulerian SPH and rigorous equivalence to finite-volume methods through particle-relaxation, kernel correction, and face-normal alignment (Wang et al., 2023).

The core principle remains the consistent, physically-motivated replacement of local derivatives by kernel-weighted, adjoint-structured, and, where necessary, corrected or generalized integral operators. These form a rigorously defined operator calculus for SPH and related meshless simulation techniques.