Add correction force for structural dynamics

docs/src/index.md

@@ -62,6 +62,14 @@ Modules = [Pixie]
 Pages = [joinpath("sph", "fluid", "viscosity.jl")]
 ```
 
+### Structural dynamics
+
+#### File penalty_force.jl
+```@autodocs
+Modules = [Pixie]
+Pages = [joinpath("sph", "solid", "penalty_force.jl")]
+```
+
 ## Util
 
 ```@autodocs

examples/solid/oscillating_beam_2d.jl

@@ -86,7 +86,8 @@ particle_container = SolidParticleContainer(particle_coordinates, particle_veloc
                                             smoothing_kernel, smoothing_length,
                                             E, nu,
                                             n_fixed_particles=n_particles_fixed,
-                                            acceleration=(0.0, -2.0))
+                                            acceleration=(0.0, -2.0),
+                                            penalty_force=PenaltyForceGanzenmueller(alpha=0.001))
 
 semi = Semidiscretization(particle_container, neighborhood_search=SpatialHashingSearch)
 tspan = (0.0, 5.0)

src/Pixie.jl

@@ -15,6 +15,7 @@ using WriteVTK: vtk_grid, MeshCell, VTKCellTypes
 include("util.jl")
 include("sph/neighborhood_search.jl") # semidiscretization/semidiscretization.jl depends on this
 include("sph/fluid/density_calculators.jl") # containers/containers.jl depends on this
+include("sph/solid/penalty_force.jl") # containers/containers.jl depends on this
 include("containers/container.jl")
 include("semidiscretization/semidiscretization.jl")
 include("sph/sph.jl")
@@ -25,6 +26,7 @@ include("setups/rectangular_tank.jl")
 
 export Semidiscretization, FluidParticleContainer, SolidParticleContainer, semidiscretize, AliveCallback, SolutionSavingCallback
 export ContinuityDensity, SummationDensity
+export PenaltyForceGanzenmueller
 export SchoenbergCubicSplineKernel, SchoenbergQuarticSplineKernel, SchoenbergQuinticSplineKernel
 export StateEquationIdealGas, StateEquationCole
 export ArtificialViscosityMonaghan

src/containers/solid_container.jl

@@ -24,7 +24,7 @@ The discretized version of this equation is given by (O’Connor & Rogers 2021):
 ```math
 \frac{\mathrm{d}\bm{v}_a}{\mathrm{d}t} = \sum_b m_{0b}
     \left( \frac{\bm{P}_a \bm{L}_{0a}}{\rho_{0a}^2} + \frac{\bm{P}_b \bm{L}_{0b}}{\rho_{0b}^2} \right)
-    \nabla_{0a} W(\bm{X}_{ab}) + \bm{g},
+    \nabla_{0a} W(\bm{X}_{ab}) + \frac{\bm{f}_a^{HG}}{m_{0a}} + \bm{g},
 ```
 with
 ```math
@@ -63,6 +63,8 @@ and
 ```
 are the Lamé coefficients, where $E$ is the Young's modulus and $\nu$ is the Poisson ratio.
 
+For the penalty force ``\bm{f}_a^{HG}`` see [`PenaltyForceGanzenmueller`](@ref)
+
 References:
 - Joseph O’Connor, Benedict D. Rogers.
   "A fluid–structure interaction model for free-surface flows and flexible structures using
@@ -74,7 +76,7 @@ References:
   In: International Journal for Numerical Methods in Engineering 48 (2000), pages 1359–1400.
   [doi: 10.1002/1097-0207](https://doi.org/10.1002/1097-0207)
 """
-struct SolidParticleContainer{NDIMS, ELTYPE<:Real, DC, K, C} <: ParticleContainer{NDIMS}
+struct SolidParticleContainer{NDIMS, ELTYPE<:Real, DC, K, PF, C} <: ParticleContainer{NDIMS}
     initial_coordinates ::Array{ELTYPE, 2} # [dimension, particle]
     current_coordinates ::Array{ELTYPE, 2} # [dimension, particle]
     initial_velocity    ::Array{ELTYPE, 2} # [dimension, particle]
@@ -89,6 +91,7 @@ struct SolidParticleContainer{NDIMS, ELTYPE<:Real, DC, K, C} <: ParticleContaine
     smoothing_kernel    ::K
     smoothing_length    ::ELTYPE
     acceleration        ::SVector{NDIMS, ELTYPE}
+    penalty_force       ::PF
     cache               ::C
 
     function SolidParticleContainer(particle_coordinates, particle_velocities,
@@ -97,7 +100,8 @@ struct SolidParticleContainer{NDIMS, ELTYPE<:Real, DC, K, C} <: ParticleContaine
                                     smoothing_kernel, smoothing_length,
                                     young_modulus, poisson_ratio;
                                     n_fixed_particles=0,
-                                    acceleration=ntuple(_ -> 0.0, size(particle_coordinates, 1)))
+                                    acceleration=ntuple(_ -> 0.0, size(particle_coordinates, 1)),
+                                    penalty_force=nothing)
         NDIMS = size(particle_coordinates, 1)
         ELTYPE = eltype(particle_masses)
         nparticles = length(particle_masses)
@@ -115,18 +119,26 @@ struct SolidParticleContainer{NDIMS, ELTYPE<:Real, DC, K, C} <: ParticleContaine
         lame_mu = 0.5 * young_modulus / (1 + poisson_ratio)
 
         # cache = create_cache(hydrodynamic_density_calculator, ELTYPE, nparticles)
-        cache = (; )
+        cache = (; create_cache(penalty_force, young_modulus, nparticles, NDIMS, ELTYPE)...)
 
         return new{NDIMS, ELTYPE, typeof(hydrodynamic_density_calculator),
-                   typeof(smoothing_kernel), typeof(cache)}(
+                   typeof(smoothing_kernel), typeof(penalty_force), typeof(cache)}(
             particle_coordinates, current_coordinates, particle_velocities, particle_masses,
             correction_matrix, pk1_corrected, particle_material_densities,
             n_moving_particles, lame_lambda, lame_mu,
             hydrodynamic_density_calculator, smoothing_kernel, smoothing_length,
-            acceleration_, cache)
+            acceleration_, penalty_force, cache)
     end
 end
 
+function create_cache(::Nothing, young_modulus, nparticles, NDIMS, ELTYPE)
+    return (; )
+end
+
+function create_cache(::PenaltyForceGanzenmueller, young_modulus, nparticles, NDIMS, ELTYPE)
+    deformation_grad    = Array{ELTYPE, 3}(undef, NDIMS, NDIMS, nparticles)
+    return (; deformation_grad, young_modulus)
+end
 
 @inline n_moving_particles(container::SolidParticleContainer) = container.n_moving_particles
 
@@ -138,6 +150,7 @@ end
 
 
 @inline get_correction_matrix(particle, container) = extract_smatrix(container.correction_matrix, particle, container)
+@inline get_deformation_gradient(particle, container) = extract_smatrix(container.cache.deformation_grad, particle, container)
 @inline get_pk1_corrected(particle, container) = extract_smatrix(container.pk1_corrected, particle, container)
 
 @inline function extract_smatrix(array, particle, container)
@@ -216,6 +229,11 @@ end
 
 
 @inline function compute_pk1_corrected(neighborhood_search, container)
+    @unpack penalty_force = container
+    compute_pk1_corrected(penalty_force, neighborhood_search, container)
+end
+
+@inline function compute_pk1_corrected(::Nothing, neighborhood_search, container)
     @unpack pk1_corrected = container
 
     @threaded for particle in eachparticle(container)
@@ -228,6 +246,27 @@ end
     end
 end
 
+@inline function compute_pk1_corrected(::PenaltyForceGanzenmueller, neighborhood_search, container)
+    @unpack pk1_corrected, cache = container
+    @unpack deformation_grad = cache
+
+    @threaded for particle in eachparticle(container)
+        J_particle = deformation_gradient(particle, neighborhood_search, container)
+
+        # store deformation gradient
+        for j in 1:ndims(container), i in 1:ndims(container)
+            deformation_grad[i, j, particle] = J_particle[i, j]
+        end
+
+        pk1_particle = pk1_stress_tensor(J_particle, container)
+        pk1_particle_corrected = pk1_particle * get_correction_matrix(particle, container)
+
+        for j in 1:ndims(container), i in 1:ndims(container)
+            pk1_corrected[i, j, particle] = pk1_particle_corrected[i, j]
+        end
+    end
+end
+
 
 # First Piola-Kirchhoff stress tensor
 function pk1_stress_tensor(particle, neighborhood_search, container)
@@ -238,13 +277,18 @@ function pk1_stress_tensor(particle, neighborhood_search, container)
     return J * S
 end
 
+function pk1_stress_tensor(J, container)
+    S = pk2_stress_tensor(J, container)
+
+    return J * S
+end
 
 # We cannot use a variable for the number of dimensions here, it has to be hardcoded
 @inline function deformation_gradient(particle, neighborhood_search, container::SolidParticleContainer{2})
     return @SMatrix [deformation_gradient(i, j, particle, neighborhood_search, container) for i in 1:2, j in 1:2]
 end
 
-@inline function deformation_gradient(particle, container::SolidParticleContainer{3})
+@inline function deformation_gradient(particle, neighborhood_search, container::SolidParticleContainer{3})
     return @SMatrix [deformation_gradient(i, j, particle, neighborhood_search, container) for i in 1:3, j in 1:3]
 end
 
@@ -288,6 +332,44 @@ end
 end
 
 
+# TODO
+@inline function calc_penalty_force!(du, particle, neighbor, initial_pos_diff,
+                                     initial_distance, container, penalty_force::PenaltyForceGanzenmueller)
+    @unpack smoothing_kernel, smoothing_length, mass,
+        material_density, current_coordinates, cache = container
+    @unpack young_modulus = cache
+
+    current_pos_diff = get_particle_coords(particle, current_coordinates, container) -
+        get_particle_coords(neighbor, current_coordinates, container)
+    current_distance = norm(current_pos_diff)
+
+    volume_particle = mass[particle] / material_density[particle]
+    volume_neighbor = mass[neighbor] / material_density[neighbor]
+
+    kernel_ = kernel(smoothing_kernel, initial_distance, smoothing_length)
+
+    J_a = get_deformation_gradient(particle, container)
+    J_b = get_deformation_gradient(neighbor, container)
+
+    eps_sum = (J_a + J_b) * initial_pos_diff - 2 * current_pos_diff
+    delta_sum = dot(eps_sum, current_pos_diff) / current_distance
+
+    dv = 0.5 * penalty_force.alpha * volume_particle * volume_neighbor *
+        kernel_ / initial_distance^2 * young_modulus * delta_sum *
+        current_pos_diff / current_distance
+
+    for i in 1:ndims(container)
+        du[ndims(container) + i, particle] += dv[i] / mass[particle]
+    end
+
+    return du
+end
+
+@inline function calc_penalty_force!(du, particle, neighbor, initial_pos_diff,
+                                     initial_distance, container, ::Nothing)
+    return du
+end
+
 function write_variables!(u0, container::SolidParticleContainer)
     @unpack initial_coordinates, initial_velocity = container
 

src/sph/solid/penalty_force.jl

@@ -0,0 +1,43 @@
+@doc raw"""
+    PenaltyForceGanzenmueller(; alpha=0.1)
+
+In FEM there are hour-glassing correction techniques for underintegrated (reduced number of integration points)
+finite elements. This is, the element deforms without an associated increase of energy.
+There is an analogy to SPH in the sense that the SPH approximation of the deformation gradient ``\bm{J}``
+is not unique with respect to the positions of the integration points.
+In order to circumvent this, (Georg C. Ganzenmüller 2015) introduced a so-called hourglass correction force or penalty force ``f^{HG}``,
+which is given by
+```math
+\bm{f}_a^{HG} = \frac{1}{2} \alpha \sum_b \frac{V_{0a} V_{0b} W_{0ab}}{X_{ab}^2}
+                \left( E \delta_{ab}^a + E \delta_{ba}^b \right) \frac{\bm{x}_{ab}}{x_{ab}}
+```
+where ``V_{0a}`` and ``V_{0b}`` is the unit voume of particle ``a`` and ``b`` respectively and
+``W_{0ab}`` is a smoothing kernel. The subindex ``0`` indicates the reference configuration.
+``X_{ab}`` is the distance of particle ``a`` and ``b`` in the reference configuration and ``x_{ab}`` is
+the distance in the current configuration.
+
+This correction force is based on the potential energy density of a Hookean material.
+Thus, ``E`` is the Young's modulus and ``\alpha`` is a dimensionless coefficient that controls
+the amplitude of hourglass correction.
+The separation vector $\delta_{ab}^a$ indicates the change of distance which the particle separation should attain
+in order to minimize the error and is given by
+```math
+\delta_{ab}^a = \frac{\bm{\epsilon}_{ab}^a\bm{x_{ab}}}{x_{ab}}
+```
+where the error vector is defined as
+```math
+\bm{\epsilon}_{ab}^a = \bm{J}_a \bm{X}_{ab} - \bm{x}_{ab} .
+```
+
+References:
+- Georg C. Ganzenmüller.
+  "An hourglass control algorithm for Lagrangian Smooth Particle Hydrodynamics".
+  In: Computer Methods in Applied Mechanics and Engineering 286 (2015).
+  [doi: 10.1016/j.cma.2014.12.005](https://doi.org/10.1016/j.cma.2014.12.005)
+"""
+struct PenaltyForceGanzenmueller{ELTYPE}
+    alpha   :: ELTYPE
+    function PenaltyForceGanzenmueller(; alpha=1.0)
+        new{typeof(alpha)}(alpha)
+    end
+end

src/sph/solid/rhs.jl

@@ -2,7 +2,7 @@
 function interact!(du, u_particle_container, u_neighbor_container, neighborhood_search,
                    particle_container::SolidParticleContainer,
                    neighbor_container::SolidParticleContainer)
-    @unpack smoothing_kernel, smoothing_length = particle_container
+    @unpack smoothing_kernel, smoothing_length, penalty_force = particle_container
 
     # Different solids do not interact with each other (yet)
     if particle_container !== neighbor_container
@@ -23,6 +23,8 @@ function interact!(du, u_particle_container, u_neighbor_container, neighborhood_
             if sqrt(eps()) < distance <= compact_support(smoothing_kernel, smoothing_length)
                 calc_dv!(du, particle, neighbor, pos_diff, distance,
                          particle_container, neighbor_container)
+                calc_penalty_force!(du, particle, neighbor, pos_diff,
+                                    distance, particle_container, penalty_force)
             end
         end
     end