Implement bond-associated quadrature Peridynamics model for TLSPH

docs/src/general/smoothing_kernels.md

@@ -13,6 +13,7 @@ The following smoothing kernels are currently available:
 | [`Poly6Kernel`](@ref)                     | $[0, 1h]$         | $1.5$ to $2.5$        | Academic                | +         |
 | [`SpikyKernel`](@ref)                     | $[0, 1h]$         | $1.5$ to $3.0$        | Academic                | +         |
 | [`LaguerreGaussKernel`](@ref)             | $[0, 2h]$         | $1.3$ to $1.5$        | General                 | ++++      |
+| [`ParabolicKernel`](@ref)                 | $[0, 1h]$         | $2.1$ to $3.0$        | TLSPH only              | -         |
 
 Any Kernel with a stability rating of more than '+++' doesn't suffer from pairing-instability.
 
@@ -40,6 +41,7 @@ other_kernels = [
     ("Gaussian",       GaussianKernel{2}()),
     ("Poly6",          Poly6Kernel{2}()),
     ("Laguerre-Gauss", LaguerreGaussKernel{2}()),
+    ("Parabolic",      ParabolicKernel{2}()),
 ]
 
 spiky_kernel_group = [

docs/src/refs.bib

@@ -342,6 +342,19 @@ @Article{Ganzenmueller2015
   publisher = {Elsevier {BV}},
 }
 
+@Article{Ganzenmueller2015b,
+  author    = {Ganzenm{\"{u}}ller, Georg C. and Hiermaier, Stefan and May, Michael},
+  title     = {On the similarity of meshless discretizations of {Peridynamics} and {Smooth-Particle} {Hydrodynamics}},
+  journal   = {Computers \& Structures},
+  year      = {2015},
+  volume    = {150},
+  pages     = {71--78},
+  month     = apr,
+  issn      = {0045-7949},
+  doi       = {10.1016/j.compstruc.2014.12.011},
+  publisher = {Elsevier {BV}},
+}
+
 @Book{Gasser2021,
   author    = {Gasser, T. Christian},
   title     = {Vascular Biomechanics: Concepts, Models, and Applications},

docs/src/systems/total_lagrangian_sph.md

@@ -57,6 +57,57 @@ are the Lamé coefficients, where $E$ is the Young's modulus and $\nu$ is the Po
 
 The term $\bm{f}_a^{PF}$ is an optional penalty force. See e.g. [`PenaltyForceGanzenmueller`](@ref).
 
+### Bond-Associated Quadrature
+
+The traditional formulation above is selected with `model=StandardTLSPHModel()`.
+With `model=BondAssociatedTLSPHModel()`, the corrected gradient is interpreted as the
+`:C1` reproducing-kernel gradient weight used by the bond-associated correspondence
+formulation in Peridynamics.jl. For
+$\bm{\xi}_{ab} = \bm{X}_b - \bm{X}_a$ and particle volume
+$V_b = m_{0b}/\rho_{0b}$, define
+```math
+\omega_{ab} =
+-\frac{\nabla_{0a} W(\bm{X}_{ab}) \cdot \bm{X}_{ab}}
+      {\lVert\bm{X}_{ab}\rVert^2}, \qquad
+w_a = \sum_b \omega_{ab} V_b,
+```
+and
+```math
+\bm{\Phi}_{ab} = V_b \bm{L}_{0a} \nabla_{0a} W(\bm{X}_{ab}).
+```
+These are algebraically the `:C1` moment matrix and gradient weights.
+
+The deformation gradient associated with the center of bond $ab$ is
+```math
+\bm{F}_{ab}
+= \frac{\bm{F}_a + \bm{F}_b}{2}
++ \left(\bm{x}_b - \bm{x}_a
+  - \frac{\bm{F}_a + \bm{F}_b}{2}\bm{\xi}_{ab}\right)
+  \frac{\bm{\xi}_{ab}^T}{\lVert\bm{\xi}_{ab}\rVert^2}.
+```
+Let $\bm{P}_{ab}$ be the PK1 stress evaluated from $\bm{F}_{ab}$ with the material
+parameters of particle $a$. The transverse stress integral is
+```math
+\bm{Q}_a = \sum_b \omega_{ab}
+\left(\frac{1}{2w_a} + \frac{1}{2w_b}\right) V_b
+\bm{P}_{ab}
+\left(\bm{I} - \frac{\bm{\xi}_{ab}\bm{\xi}_{ab}^T}
+{\lVert\bm{\xi}_{ab}\rVert^2}\right).
+```
+The directed force state and acceleration are then
+```math
+\bm{t}_{ab}
+= \frac{\omega_{ab}}{w_a\lVert\bm{\xi}_{ab}\rVert^2}
+  \bm{P}_{ab}\bm{\xi}_{ab}
++ \frac{\bm{Q}_a\bm{\Phi}_{ab}}{V_b},
+\qquad
+\frac{\mathrm{d}\bm{v}_a}{\mathrm{d}t}
+= \frac{1}{\rho_{0a}}\sum_b V_b(\bm{t}_{ab} - \bm{t}_{ba}).
+```
+This is the bond-associated quadrature update implemented in
+[`Peridynamics.jl`](https://github.com/kaipartmann/Peridynamics.jl/blob/main/src/physics/rk_correspondence.jl),
+adapted to the existing TLSPH kernels and two- or three-dimensional particle systems.
+
 ```@autodocs
 Modules = [TrixiParticles]
 Pages = [joinpath("schemes", "structure", "total_lagrangian_sph", "system.jl")]

examples/structure/oscillating_beam_2d.jl

@@ -67,6 +67,7 @@ structure_system = TotalLagrangianSPHSystem(structure; smoothing_kernel, smoothi
                                             acceleration=(0.0, -gravity),
                                             penalty_force=nothing, viscosity=nothing,
                                             clamped_particles_motion=nothing,
+                                            model=StandardTLSPHModel(),
                                             self_interaction_nhs=:default)
 
 # ==========================================================================================

src/TrixiParticles.jl

@@ -67,6 +67,7 @@ include("visualization/recipes_plots.jl")
 export Semidiscretization, semidiscretize, restart_with!
 export InitialCondition, apply_angular_velocity
 export WeaklyCompressibleSPHSystem, EntropicallyDampedSPHSystem, TotalLagrangianSPHSystem,
+       StandardTLSPHModel, BondAssociatedTLSPHModel,
        RigidBodySystem, WallBoundarySystem, DEMSystem, BoundaryDEMSystem,
        OpenBoundarySystem,
        ImplicitIncompressibleSPHSystem
@@ -82,7 +83,7 @@ export PenaltyForceGanzenmueller, TransportVelocityAdami, ParticleShiftingTechni
        ContinuityEquationTermSun2019, MomentumEquationTermSun2019, VelocityAveraging
 export SchoenbergCubicSplineKernel, SchoenbergQuarticSplineKernel,
        SchoenbergQuinticSplineKernel, GaussianKernel, WendlandC2Kernel, WendlandC4Kernel,
-       WendlandC6Kernel, SpikyKernel, Poly6Kernel, LaguerreGaussKernel
+       WendlandC6Kernel, SpikyKernel, Poly6Kernel, LaguerreGaussKernel, ParabolicKernel
 export StateEquationCole, StateEquationIdealGas, StateEquationAdaptiveCole
 export ArtificialViscosityMonaghan, ViscosityAdami, ViscosityMorris, ViscosityAdamiSGS,
        ViscosityMorrisSGS, ViscosityCarreauYasuda

src/general/smoothing_kernels.jl

@@ -820,3 +820,68 @@ end
     # C' = 1 / (4 * pi * h^3 * C) = C_ / (h^3)
     return C_ * h_inv^3
 end
+
+@doc raw"""
+    ParabolicKernel{NDIMS}()
+
+Parabolic smoothing kernel given by
+```math
+    W(r, h) = \frac{1}{h^d} w(r/h)
+```
+with
+```math
+w(q) = \sigma \begin{cases}
+    1 - q^2   & \text{if } 0 \leq q < 1, \\
+    0         & \text{if } q \geq 1,
+\end{cases}
+```
+where ``d`` is the number of dimensions and ``\sigma`` is a normalization factor.
+The normalization factor ``\sigma`` is ``3 / 4`` in one dimension, ``2 / \pi``
+in two dimensions, and ``15 / (8 \pi)`` in three dimensions.
+
+This kernel function has a compact support of ``[0, h]``.
+
+The parabolic kernel is intended for Total Lagrangian SPH (TLSPH), where the
+gradient is evaluated with respect to initial coordinates. With
+``X_{ij} = X_i - X_j``, this kernel yields
+```math
+    \nabla W_{ij} = -c X_{ij},
+```
+where ``c`` is constant for fixed smoothing length. This is the kernel that makes
+the basic Peridynamics discretization equivalent to TLSPH, as shown by
+[Ganzenmüller et al. (2015)](@cite Ganzenmueller2015b).
+
+!!! warning
+    This kernel is not useful for fluid SPH. Its simple gradient causes severe
+    particle clustering in fluid simulations. Use it only for TLSPH applications.
+
+For general information and usage see [Smoothing Kernels](@ref smoothing_kernel).
+"""
+struct ParabolicKernel{NDIMS} <: AbstractSmoothingKernel{NDIMS} end
+
+@inline function kernel_unsafe(kernel::ParabolicKernel, r::Real, h)
+    h_inv = 1 / h
+    q = r * h_inv
+
+    return normalization_factor(kernel, h_inv) * (1 - q^2)
+end
+
+@inline function kernel_deriv_div_r_unsafe(kernel::ParabolicKernel, r::Real, h)
+    h_inv = 1 / h
+
+    return -2 * normalization_factor(kernel, h_inv) * h_inv^2
+end
+
+@inline compact_support(::ParabolicKernel, h) = h
+
+@inline function normalization_factor(::ParabolicKernel{1}, h_inv)
+    return oftype(h_inv, 3 / 4) * h_inv
+end
+
+@inline function normalization_factor(::ParabolicKernel{2}, h_inv)
+    return oftype(h_inv, 2 / pi) * h_inv^2
+end
+
+@inline function normalization_factor(::ParabolicKernel{3}, h_inv)
+    return oftype(h_inv, 15 / (8 * pi)) * h_inv^3
+end

src/schemes/structure/total_lagrangian_sph/rhs.jl

@@ -18,6 +18,13 @@ end
 # Function barrier without dispatch for unit testing
 @inline function interact_structure_structure!(dv, v_system, system, semi;
                                                eachparticle=each_integrated_particle(system))
+    interact_structure_structure!(dv, v_system, system, tlsph_model(system), semi;
+                                  eachparticle)
+end
+
+@inline function interact_structure_structure!(dv, v_system, system,
+                                               ::StandardTLSPHModel, semi;
+                                               eachparticle=each_integrated_particle(system))
     (; penalty_force) = system
 
     # Everything here is done in the initial coordinates
@@ -101,6 +108,95 @@ end
     return dv
 end
 
+@inline function interact_structure_structure!(dv, v_system, system,
+                                               ::BondAssociatedTLSPHModel, semi;
+                                               eachparticle=each_integrated_particle(system))
+    (; penalty_force, mass, material_density) = system
+    (; weighted_volume) = system.cache
+
+    system_coords = initial_coordinates(system)
+    neighborhood_search = get_neighborhood_search(system, semi)
+    backend = semi.parallelization_backend
+    h = initial_smoothing_length(system)
+    almostzero = sqrt(eps(h^2))
+
+    @threaded semi for particle in eachparticle
+        m_a = @inbounds mass[particle]
+        rho_a = @inbounds material_density[particle]
+        volume_a = div_fast(m_a, rho_a)
+        current_coords_a = @inbounds current_coords(system, particle)
+        F_a = @inbounds deformation_gradient(system, particle)
+        stress_integral_a = @inbounds bond_associated_stress_integral(system, particle)
+        weighted_volume_a = @inbounds weighted_volume[particle]
+
+        dv_particle = Ref(zero(current_coords_a))
+
+        @inbounds foreach_neighbor(system_coords, system_coords,
+                                   neighborhood_search, backend,
+                                   particle) do particle, neighbor,
+                                                initial_pos_diff, initial_distance
+            initial_distance < almostzero && return
+
+            rho_b = material_density[neighbor]
+            m_b = mass[neighbor]
+            volume_b = div_fast(m_b, rho_b)
+            current_coords_b = current_coords(system, neighbor)
+            F_b = deformation_gradient(system, neighbor)
+
+            current_pos_diff_ = current_coords_a - current_coords_b
+            current_pos_diff = convert.(eltype(system), current_pos_diff_)
+            current_distance = norm(current_pos_diff)
+            initial_bond = -initial_pos_diff
+            current_bond = -current_pos_diff
+
+            F_ab = bond_deformation_gradient(F_a, F_b, initial_bond, current_bond,
+                                             initial_distance)
+            P_ab = pk1_stress_tensor(F_ab, system, particle)
+            P_ba = pk1_stress_tensor(F_ab, system, neighbor)
+
+            grad_kernel_a = smoothing_kernel_grad_unsafe(system, initial_pos_diff,
+                                                         initial_distance, particle)
+            weight_a = bond_weight(grad_kernel_a, initial_pos_diff, initial_distance)
+            gradient_weight_a = volume_b * correction_matrix(system, particle) *
+                                grad_kernel_a
+
+            neighbor_pos_diff = -initial_pos_diff
+            grad_kernel_b = smoothing_kernel_grad_unsafe(system, neighbor_pos_diff,
+                                                         initial_distance, neighbor)
+            weight_b = bond_weight(grad_kernel_b, neighbor_pos_diff, initial_distance)
+            gradient_weight_b = volume_a * correction_matrix(system, neighbor) *
+                                grad_kernel_b
+
+            stress_integral_b = bond_associated_stress_integral(system, neighbor)
+            weighted_volume_b = weighted_volume[neighbor]
+
+            force_state_ab = weight_a / (weighted_volume_a * initial_distance^2) *
+                             (P_ab * initial_bond) +
+                             stress_integral_a * gradient_weight_a / volume_b
+            force_state_ba = weight_b / (weighted_volume_b * initial_distance^2) *
+                             (P_ba * -initial_bond) +
+                             stress_integral_b * gradient_weight_b / volume_a
+
+            dv_particle[] += volume_b / rho_a * (force_state_ab - force_state_ba)
+
+            @inbounds dv_penalty_force!(dv_particle, penalty_force, particle, neighbor,
+                                        initial_pos_diff, initial_distance,
+                                        current_pos_diff, current_distance,
+                                        system, m_a, m_b, rho_a, rho_b, F_a, F_b)
+
+            @inbounds dv_viscosity_tlsph!(dv_particle, system, v_system, particle, neighbor,
+                                          current_pos_diff, current_distance,
+                                          m_a, m_b, rho_a, rho_b, F_a, grad_kernel_a)
+        end
+
+        for i in 1:ndims(system)
+            @inbounds dv[i, particle] += dv_particle[][i]
+        end
+    end
+
+    return dv
+end
+
 # Structure-fluid interaction
 function interact!(dv, v_particle_system, u_particle_system,
                    v_neighbor_system, u_neighbor_system,