Merge branch 'main' into dev

Author: LasNikas

docs/src/refs.bib

@@ -108,7 +108,7 @@ @article{Balsara1995
 pages = {357-372},
 year = {1995},
 issn = {0021-9991},
-doi = {https://doi.org/10.1016/S0021-9991(95)90221-X},
+doi = {10.1016/S0021-9991(95)90221-X},
 url = {https://www.sciencedirect.com/science/article/pii/S002199919590221X},
 author = {Dinshaw S. Balsara},
 }
@@ -406,6 +406,18 @@ @Article{Li1996
   publisher = {Elsevier BV},
 }
 
+@article{Lin2015,
+	author = {Lin, Jun and Naceur, Hakim and Coutellier, Daniel and Laksimi, Abdel},
+	journal = {Engineering Computations},
+	title = {Geometrically nonlinear analysis of two-dimensional structures using an improved smoothed particle hydrodynamics method},
+	year = {2015},
+	issn = {0264-4401},
+  volume = {32},
+	number = {3},
+  pages = {779--805},
+	doi = {10.1108/EC-12-2013-0306},
+}
+
 @Article{Liu2006,
   author    = {Liu, M.B. and Liu, G.R.},
   journal   = {Applied Numerical Mathematics},
@@ -544,7 +556,7 @@ @article{Morris2000
   number = {3},
   pages = {333-353},
   keywords = {interfacial flow, meshless methods, surface tension},
-  doi = {https://doi.org/10.1002/1097-0363(20000615)33:3<333::AID-FLD11>3.0.CO;2-7},
+  doi = {10.1002/1097-0363(20000615)33:3<333::AID-FLD11>3.0.CO;2-7},
   year = {2000}
 }
 
@@ -735,8 +747,6 @@ @Article{Wendland1995
   publisher = {Springer Science and Business Media LLC},
 }
 
-@Comment{jabref-meta: databaseType:bibtex;}
-
 @book{Lange2005,
   author = {{Lange's Handbook of Chemistry}},
   title = {Lange's Handbook of Chemistry},

docs/src/systems/fluid.md

@@ -55,7 +55,7 @@ by Balsara ([Balsara1995](@cite)) or Morris ([Morris1997](@cite)).
 
 ##### Mathematical Formulation
 
-The force exerted by particle `b` on particle `a` due to artificial viscosity is given by:
+The force exerted by particle ``b`` on particle ``a`` due to artificial viscosity is given by:
 
 ```math
 F_{ab}^{\text{AV}} = - m_a m_b \Pi_{ab} \nabla W_{ab}

docs/src/systems/total_lagrangian_sph.md

@@ -102,3 +102,30 @@ where the error vector is defined as
 Modules = [TrixiParticles]
 Pages = [joinpath("schemes", "solid", "total_lagrangian_sph", "penalty_force.jl")]
 ```
+
+## Viscosity
+
+Another technique that is used to correct the hourglass instability is artificial viscosity.
+Hereby, a viscosity term designed for fluids (see [Viscosity](@ref viscosity_sph)) is applied.
+First, the force ``f_{ab}^{\text{fluid}}`` exerted by particle ``b`` on particle ``a``
+due to artificial viscosity is computed as if both particles were fluid particles
+(see [Viscosity](@ref viscosity_sph) for the relevant equations).
+Then, according to [Lin et al. (2015)](@cite Lin2015), this force can be applied to TLSPH
+with the following conversion:
+```math
+f_{ab}^{\text{AV}} = \det(F_a) F_a^{-1} f_{ab}^{\text{fluid}},
+```
+where ``F_a`` is the deformation gradient at particle ``a``.
+
+We found that artificial viscosity is not effective at correcting the incorrect
+particle positions due to hourglass modes.
+It does however prevent particles from oscillating between different incorrect positions,
+which develops into an instability if uncorrected.
+We still recommend penalty force over artificial viscosity to correct hourglass modes
+as penalty force is specifically designed to correct the incorrect particle positions.
+
+In some FSI simulations, notably when very thin structures or structures with low material
+density are present, instabilities in the fluid can be induced by the structure.
+In these cases, artificial viscosity is effective at stabilizing the fluid close to the
+structure, and we recommend using it in combination with penalty force to both
+prevent hourglass modes and stabilize the fluid close to the fluid-structure interface.

examples/solid/oscillating_beam_2d.jl

@@ -63,7 +63,7 @@ solid_system = TotalLagrangianSPHSystem(solid, smoothing_kernel, smoothing_lengt
                                         material.E, material.nu,
                                         n_fixed_particles=nparticles(fixed_particles),
                                         acceleration=(0.0, -gravity),
-                                        penalty_force=nothing)
+                                        penalty_force=nothing, viscosity=nothing)
 
 # ==========================================================================================
 # ==== Simulation

src/schemes/fluid/viscosity.jl

@@ -1,5 +1,6 @@
-
-# Unpack the neighboring systems viscosity to dispatch on the viscosity type
+# Unpack the neighboring systems viscosity to dispatch on the viscosity type.
+# This function is only necessary to allow `nothing` as viscosity.
+# Otherwise, we could just apply the viscosity as a function directly.
 @propagate_inbounds function dv_viscosity(particle_system, neighbor_system,
                                           v_particle_system, v_neighbor_system,
                                           particle, neighbor, pos_diff, distance,
@@ -98,6 +99,7 @@ end
 
     smoothing_length_particle = smoothing_length(particle_system, particle)
     smoothing_length_neighbor = smoothing_length(particle_system, neighbor)
+    smoothing_length_average = (smoothing_length_particle + smoothing_length_neighbor) / 2
 
     nu_a = kinematic_viscosity(particle_system,
                                viscosity_model(neighbor_system, particle_system),
@@ -106,7 +108,6 @@ end
                                viscosity_model(particle_system, neighbor_system),
                                smoothing_length_neighbor, sound_speed)
 
-    smoothing_length_average = (smoothing_length_particle + smoothing_length_neighbor) / 2
     pi_ab = viscosity(sound_speed, v_diff, pos_diff, distance, rho_mean, rho_a, rho_b,
                       smoothing_length_average, grad_kernel, nu_a, nu_b)
 

src/schemes/solid/total_lagrangian_sph/penalty_force.jl

@@ -14,45 +14,35 @@ struct PenaltyForceGanzenmueller{ELTYPE}
     end
 end
 
-@inline function calc_penalty_force!(dv, particle, neighbor, initial_pos_diff,
-                                     initial_distance, system, m_a, m_b, rho_a, rho_b,
-                                     penalty_force::PenaltyForceGanzenmueller)
-    (; young_modulus) = system
-
-    current_pos_diff = current_coords(system, particle) -
-                       current_coords(system, neighbor)
-    current_distance = norm(current_pos_diff)
+@inline function dv_penalty_force(penalty_force::Nothing,
+                                  particle, neighbor, initial_pos_diff, initial_distance,
+                                  current_pos_diff, current_distance,
+                                  system, m_a, m_b, rho_a, rho_b)
+    return zero(initial_pos_diff)
+end
 
-    volume_particle = m_a / rho_a
-    volume_neighbor = m_b / rho_b
+@inline function dv_penalty_force(penalty_force::PenaltyForceGanzenmueller,
+                                  particle, neighbor, initial_pos_diff, initial_distance,
+                                  current_pos_diff, current_distance,
+                                  system, m_a, m_b, rho_a, rho_b)
+    volume_a = m_a / rho_a
+    volume_b = m_b / rho_b
 
     kernel_weight = smoothing_kernel(system, initial_distance, particle)
 
-    J_a = deformation_gradient(system, particle)
-    J_b = deformation_gradient(system, neighbor)
+    F_a = deformation_gradient(system, particle)
+    F_b = deformation_gradient(system, neighbor)
 
     # Use the symmetry of epsilon to simplify computations
-    eps_sum = (J_a + J_b) * initial_pos_diff - 2 * current_pos_diff
+    eps_sum = (F_a + F_b) * initial_pos_diff - 2 * current_pos_diff
     delta_sum = dot(eps_sum, current_pos_diff) / current_distance
 
-    E = young_modulus_per_particle(young_modulus, particle)
+    E = young_modulus(system, particle)
 
-    f = (penalty_force.alpha / 2) * volume_particle * volume_neighbor *
+    f = (penalty_force.alpha / 2) * volume_a * volume_b *
         kernel_weight / initial_distance^2 * E * delta_sum * current_pos_diff /
         current_distance
 
-    @inbounds for i in 1:ndims(system)
-        # Divide force by mass to obtain acceleration
-        dv[i, particle] += f[i] / m_a
-    end
-
-    return dv
-end
-
-function young_modulus_per_particle(young_modulus::AbstractVector, particle)
-    return young_modulus[particle]
-end
-
-function young_modulus_per_particle(young_modulus, particle)
-    return young_modulus
+    # Divide force by mass to obtain acceleration
+    return f / m_a
 end

src/schemes/solid/total_lagrangian_sph/rhs.jl

@@ -3,56 +3,61 @@ function interact!(dv, v_particle_system, u_particle_system,
                    v_neighbor_system, u_neighbor_system,
                    particle_system::TotalLagrangianSPHSystem,
                    neighbor_system::TotalLagrangianSPHSystem, semi)
-    interact_solid_solid!(dv, particle_system, neighbor_system, semi)
+    # Different solids do not interact with each other (yet)
+    particle_system !== neighbor_system && return dv
+
+    interact_solid_solid!(dv, v_particle_system, particle_system, semi)
 end
 
 # Function barrier without dispatch for unit testing
-@inline function interact_solid_solid!(dv, particle_system, neighbor_system, semi)
-    (; penalty_force) = particle_system
-
-    # Different solids do not interact with each other (yet)
-    if particle_system !== neighbor_system
-        return dv
-    end
+@inline function interact_solid_solid!(dv, v_system, system, semi)
+    (; penalty_force) = system
 
     # Everything here is done in the initial coordinates
-    system_coords = initial_coordinates(particle_system)
-    neighbor_coords = initial_coordinates(neighbor_system)
+    system_coords = initial_coordinates(system)
 
     # Loop over all pairs of particles and neighbors within the kernel cutoff.
     # For solid-solid interaction, this has to happen in the initial coordinates.
-    foreach_point_neighbor(particle_system, neighbor_system, system_coords, neighbor_coords,
-                           semi;
-                           points=each_moving_particle(particle_system)) do particle,
-                                                                            neighbor,
-                                                                            initial_pos_diff,
-                                                                            initial_distance
+    foreach_point_neighbor(system, system, system_coords, system_coords, semi;
+                           points=each_moving_particle(system)) do particle, neighbor,
+                                                                   initial_pos_diff,
+                                                                   initial_distance
         # Only consider particles with a distance > 0.
         initial_distance < sqrt(eps()) && return
 
-        rho_a = particle_system.material_density[particle]
-        rho_b = neighbor_system.material_density[neighbor]
+        rho_a = @inbounds system.material_density[particle]
+        rho_b = @inbounds system.material_density[neighbor]
 
-        grad_kernel = smoothing_kernel_grad(particle_system, initial_pos_diff,
+        grad_kernel = smoothing_kernel_grad(system, initial_pos_diff,
                                             initial_distance, particle)
 
-        m_a = particle_system.mass[particle]
-        m_b = neighbor_system.mass[neighbor]
+        m_a = @inbounds system.mass[particle]
+        m_b = @inbounds system.mass[neighbor]
 
-        dv_particle = m_b *
-                      (pk1_corrected(particle_system, particle) / rho_a^2 +
-                       pk1_corrected(neighbor_system, neighbor) / rho_b^2) *
-                      grad_kernel
+        current_pos_diff = @inbounds current_coords(system, particle) -
+                                     current_coords(system, neighbor)
+        current_distance = norm(current_pos_diff)
 
-        @inbounds for i in 1:ndims(particle_system)
-            dv[i, particle] += dv_particle[i]
-        end
+        dv_stress = m_b *
+                    (pk1_corrected(system, particle) / rho_a^2 +
+                     pk1_corrected(system, neighbor) / rho_b^2) * grad_kernel
 
-        calc_penalty_force!(dv, particle, neighbor, initial_pos_diff,
-                            initial_distance, particle_system, m_a, m_b, rho_a, rho_b,
-                            penalty_force)
+        dv_penalty_force_ = dv_penalty_force(penalty_force, particle, neighbor,
+                                             initial_pos_diff, initial_distance,
+                                             current_pos_diff, current_distance,
+                                             system, m_a, m_b, rho_a, rho_b)
+
+        dv_viscosity = dv_viscosity_tlsph(system, v_system, particle, neighbor,
+                                          current_pos_diff, current_distance,
+                                          m_a, m_b, rho_a, rho_b, grad_kernel)
+
+        for i in 1:ndims(system)
+            @inbounds dv[i,
+                         particle] += dv_stress[i] + dv_penalty_force_[i] +
+                                      dv_viscosity[i]
+        end
 
-        # TODO continuity equation?
+        # TODO continuity equation for boundary model with `ContinuityDensity`?
     end
 
     return dv
@@ -90,10 +95,10 @@ function interact!(dv, v_particle_system, u_particle_system,
 
         rho_a = current_density(v_particle_system, particle_system, particle)
         rho_b = current_density(v_neighbor_system, neighbor_system, neighbor)
-        rho_mean = (rho_a + rho_b) / 2
 
         # Use kernel from the fluid system in order to get the same force here in
         # solid-fluid interaction as for fluid-solid interaction.
+        # TODO this will not use corrections if the fluid uses corrections.
         grad_kernel = smoothing_kernel_grad(neighbor_system, pos_diff, distance, particle)
 
         # In fluid-solid interaction, use the "hydrodynamic pressure" of the solid particles

src/schemes/solid/total_lagrangian_sph/system.jl

@@ -30,6 +30,8 @@ See [Total Lagrangian SPH](@ref tlsph) for more details on the method.
                     fluid-structure interaction (see [Boundary Models](@ref boundary_models)).
 - `penalty_force`:  Penalty force to ensure regular particle position under large deformations
                     (see [`PenaltyForceGanzenmueller`](@ref)).
+- `viscosity`:      Artificial viscosity model to stabilize both the TLSPH and the FSI.
+                    Currently, only [`ArtificialViscosityMonaghan`](@ref) is supported.
 - `acceleration`:   Acceleration vector for the system. (default: zero vector)
 - `source_terms`:   Additional source terms for this system. Has to be either `nothing`
                     (by default), or a function of `(coords, velocity, density, pressure)`
@@ -55,7 +57,7 @@ See [Total Lagrangian SPH](@ref tlsph) for more details on the method.
     where `beam` and `fixed_particles` are of type `InitialCondition`.
 """
 struct TotalLagrangianSPHSystem{BM, NDIMS, ELTYPE <: Real, IC, ARRAY1D, ARRAY2D, ARRAY3D,
-                                YM, PR, LL, LM, K, PF, ST} <: SolidSystem{NDIMS}
+                                YM, PR, LL, LM, K, PF, V, ST} <: SolidSystem{NDIMS}
     initial_condition   :: IC
     initial_coordinates :: ARRAY2D # Array{ELTYPE, 2}: [dimension, particle]
     current_coordinates :: ARRAY2D # Array{ELTYPE, 2}: [dimension, particle]
@@ -74,6 +76,7 @@ struct TotalLagrangianSPHSystem{BM, NDIMS, ELTYPE <: Real, IC, ARRAY1D, ARRAY2D,
     acceleration        :: SVector{NDIMS, ELTYPE}
     boundary_model      :: BM
     penalty_force       :: PF
+    viscosity           :: V
     source_terms        :: ST
     buffer              :: Nothing
 end
@@ -84,7 +87,8 @@ function TotalLagrangianSPHSystem(initial_condition,
                                   n_fixed_particles=0, boundary_model=nothing,
                                   acceleration=ntuple(_ -> 0.0,
                                                       ndims(smoothing_kernel)),
-                                  penalty_force=nothing, source_terms=nothing)
+                                  penalty_force=nothing, viscosity=nothing,
+                                  source_terms=nothing)
     NDIMS = ndims(initial_condition)
     ELTYPE = eltype(initial_condition)
     n_particles = nparticles(initial_condition)
@@ -119,7 +123,7 @@ function TotalLagrangianSPHSystem(initial_condition,
                                     n_moving_particles, young_modulus, poisson_ratio,
                                     lame_lambda, lame_mu, smoothing_kernel,
                                     smoothing_length, acceleration_, boundary_model,
-                                    penalty_force, source_terms, nothing)
+                                    penalty_force, viscosity, source_terms, nothing)
 end
 
 function Base.show(io::IO, system::TotalLagrangianSPHSystem)
@@ -130,6 +134,7 @@ function Base.show(io::IO, system::TotalLagrangianSPHSystem)
     print(io, ", ", system.acceleration)
     print(io, ", ", system.boundary_model)
     print(io, ", ", system.penalty_force)
+    print(io, ", ", system.viscosity)
     print(io, ") with ", nparticles(system), " particles")
 end
 
@@ -159,7 +164,8 @@ function Base.show(io::IO, ::MIME"text/plain", system::TotalLagrangianSPHSystem)
         summary_line(io, "smoothing kernel", system.smoothing_kernel |> typeof |> nameof)
         summary_line(io, "acceleration", system.acceleration)
         summary_line(io, "boundary model", system.boundary_model)
-        summary_line(io, "penalty force", system.penalty_force |> typeof |> nameof)
+        summary_line(io, "penalty force", system.penalty_force)
+        summary_line(io, "viscosity", system.viscosity)
         summary_footer(io)
     end
 end
@@ -230,6 +236,32 @@ end
     extract_smatrix(system.pk1_corrected, system, particle)
 end
 
+function young_modulus(system::TotalLagrangianSPHSystem, particle)
+    return young_modulus(system, system.young_modulus, particle)
+end
+
+function young_modulus(::TotalLagrangianSPHSystem, young_modulus, particle)
+    return young_modulus
+end
+
+function young_modulus(::TotalLagrangianSPHSystem,
+                       young_modulus::AbstractVector, particle)
+    return young_modulus[particle]
+end
+
+function poisson_ratio(system::TotalLagrangianSPHSystem, particle)
+    return poisson_ratio(system, system.poisson_ratio, particle)
+end
+
+function poisson_ratio(::TotalLagrangianSPHSystem, poisson_ratio, particle)
+    return poisson_ratio
+end
+
+function poisson_ratio(::TotalLagrangianSPHSystem,
+                       poisson_ratio::AbstractVector, particle)
+    return poisson_ratio[particle]
+end
+
 function initialize!(system::TotalLagrangianSPHSystem, semi)
     (; correction_matrix) = system
 
@@ -345,12 +377,6 @@ end
     return lame_lambda * tr(E) * I + 2 * lame_mu * E
 end
 
-@inline function calc_penalty_force!(dv, particle, neighbor, initial_pos_diff,
-                                     initial_distance, system, m_a, m_b, rho_a, rho_b,
-                                     ::Nothing)
-    return dv
-end
-
 function write_u0!(u0, system::TotalLagrangianSPHSystem)
     (; initial_condition) = system