diff --git a/docs/src/systems/total_lagrangian_sph.md b/docs/src/systems/total_lagrangian_sph.md
index 86167d164b..d9495c08e6 100644
--- a/docs/src/systems/total_lagrangian_sph.md
+++ b/docs/src/systems/total_lagrangian_sph.md
@@ -129,3 +129,26 @@ 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.
+
+
+## [Velocity Averaging](@id velocity_averaging)
+
+In challenging FSI cases with very stiff structures, the two techniques above might not be
+sufficient to prevent instabilities.
+High-frequency noise in the structure velocity can trigger instabilities or spurious
+pressure waves (due to aliasing) in the fluid.
+Another stabilization technique is an exponential moving average (EMA) of the structure
+velocity used **only** for the fluid-structure viscous coupling (no-slip boundary condition).
+
+The averaged velocity ``\bar v`` is updated (by the [`UpdateCallback`](@ref)
+or in every sub-step of the [`SplitIntegrationCallback`](@ref) if it is used) as
+```math
+\bar v^{n+1} = (1-\alpha)\,\bar v^n + \alpha\, v^{n+1}, \qquad
+\alpha = 1 - \exp(-\Delta t/\tau),
+```
+where ``\tau`` is the time constant.
+
+```@autodocs
+Modules = [TrixiParticles]
+Pages = [joinpath("schemes", "structure", "total_lagrangian_sph", "velocity_averaging.jl")]
+```
diff --git a/examples/fsi/dam_break_plate_2d.jl b/examples/fsi/dam_break_plate_2d.jl
index 44d86c5317..176a3efd68 100644
--- a/examples/fsi/dam_break_plate_2d.jl
+++ b/examples/fsi/dam_break_plate_2d.jl
@@ -129,7 +129,8 @@ structure_system = TotalLagrangianSPHSystem(structure,
                                             E, nu, boundary_model=boundary_model_structure,
                                             clamped_particles=1:nparticles(clamped_particles),
                                             acceleration=(0.0, -gravity),
-                                            penalty_force=PenaltyForceGanzenmueller(alpha=0.01))
+                                            penalty_force=PenaltyForceGanzenmueller(alpha=0.01),
+                                            velocity_averaging=nothing)
 
 # ==========================================================================================
 # ==== Simulation
diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl
index bc3cbca3a2..ad2c07d4c5 100644
--- a/src/TrixiParticles.jl
+++ b/src/TrixiParticles.jl
@@ -76,7 +76,7 @@ export InfoCallback, SolutionSavingCallback, DensityReinitializationCallback,
 export ContinuityDensity, SummationDensity
 export PenaltyForceGanzenmueller, TransportVelocityAdami, ParticleShiftingTechnique,
        ParticleShiftingTechniqueSun2017, ConsistentShiftingSun2019,
-       ContinuityEquationTermSun2019, MomentumEquationTermSun2019
+       ContinuityEquationTermSun2019, MomentumEquationTermSun2019, VelocityAveraging
 export SchoenbergCubicSplineKernel, SchoenbergQuarticSplineKernel,
        SchoenbergQuinticSplineKernel, GaussianKernel, WendlandC2Kernel, WendlandC4Kernel,
        WendlandC6Kernel, SpikyKernel, Poly6Kernel, LaguerreGaussKernel
diff --git a/src/callbacks/split_integration.jl b/src/callbacks/split_integration.jl
index 372ec752f0..02d1543c68 100644
--- a/src/callbacks/split_integration.jl
+++ b/src/callbacks/split_integration.jl
@@ -48,11 +48,22 @@ callback = SplitIntegrationCallback(RDPK3SpFSAL49(), abstol=1e-6, reltol=1e-4)
 │ alg: …………………………………………………………………… RDPK3SpFSAL49                                                    │
 │ abstol: …………………………………………………………… 1.0e-6                                                           │
 │ reltol: …………………………………………………………… 0.0001                                                           │
+│ callback: ……………………………………………………… UpdateAveragedVelocityCallback                                   │
 └──────────────────────────────────────────────────────────────────────────────────────────────────┘
 ```
 """
 function SplitIntegrationCallback(alg; kwargs...)
-    split_integration_callback = SplitIntegrationCallback(nothing, alg, kwargs)
+    # Add lightweight callback to (potentially) update the averaged velocity
+    # during the split integration.
+    if haskey(kwargs, :callback)
+        # Note that `CallbackSet`s can be nested
+        kwargs = (; kwargs...,
+                  callback=CallbackSet(values(kwargs).callback,
+                                       UpdateAveragedVelocityCallback()))
+    else
+        kwargs = (; kwargs..., callback=UpdateAveragedVelocityCallback())
+    end
+    split_integration_callback = SplitIntegrationCallback(nothing, alg, pairs(kwargs))
 
     # The first one is the `condition`, the second the `affect!`
     return DiscreteCallback(split_integration_callback, split_integration_callback,
@@ -356,3 +367,53 @@ function Base.show(io::IO, ::MIME"text/plain",
         summary_box(io, "SplitIntegrationCallback", setup)
     end
 end
+
+# === Non-public callback for updating the averaged velocity ===
+# When no split integration is used, this is done from the `UpdateCallback`.
+# With split integration, we use this lightweight callback to avoid updating the systems.
+function UpdateAveragedVelocityCallback()
+    # The first one is the `condition`, the second the `affect!`
+    return DiscreteCallback(update_averaged_velocity_callback!,
+                            update_averaged_velocity_callback!,
+                            initialize=(initialize_averaged_velocity_callback!),
+                            save_positions=(false, false))
+end
+
+# `initialize`
+function initialize_averaged_velocity_callback!(cb, vu_ode, t, integrator)
+    v_ode, u_ode = vu_ode.x
+    semi = integrator.p.semi_split
+
+    foreach_system(semi) do system
+        initialize_averaged_velocity!(system, v_ode, semi, t)
+    end
+
+    return cb
+end
+
+# `condition`
+function update_averaged_velocity_callback!(u, t, integrator)
+    return true
+end
+
+# `affect!`
+function update_averaged_velocity_callback!(integrator)
+    t_new = integrator.t
+    semi = integrator.p.semi_split
+    v_ode, u_ode = integrator.u.x
+
+    foreach_system(semi) do system
+        compute_averaged_velocity!(system, v_ode, semi, t_new)
+    end
+
+    # Tell OrdinaryDiffEq that `integrator.u` has not been modified
+    u_modified!(integrator, false)
+
+    return integrator
+end
+
+function Base.show(io::IO,
+                   cb::DiscreteCallback{typeof(update_averaged_velocity_callback!)})
+    @nospecialize cb # reduce precompilation time
+    print(io, "UpdateAveragedVelocityCallback")
+end
diff --git a/src/callbacks/update.jl b/src/callbacks/update.jl
index ab76dfacc2..b5e8b66cdb 100644
--- a/src/callbacks/update.jl
+++ b/src/callbacks/update.jl
@@ -52,12 +52,17 @@ function initial_update!(cb, u, t, integrator)
     initial_update!(cb.affect!, u, t, integrator)
 end
 
-function initial_update!(cb::UpdateCallback, u, t, integrator)
+function initial_update!(cb::UpdateCallback, vu_ode, t, integrator)
+    v_ode, u_ode = vu_ode.x
     semi = integrator.p
 
     # Tell the semidiscretization that the `UpdateCallback` is used
     semi.update_callback_used[] = true
 
+    foreach_system(semi) do system
+        initialize_averaged_velocity!(system, v_ode, semi, t)
+    end
+
     return cb(integrator)
 end
 
@@ -105,6 +110,10 @@ function (update_callback!::UpdateCallback)(integrator)
             particle_shifting_from_callback!(u_ode, shifting_technique(system), system,
                                              v_ode, semi, integrator)
         end
+
+        foreach_system(semi) do system
+            compute_averaged_velocity!(system, v_ode, semi, t)
+        end
     end
 
     return integrator
diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl
index 38531b37ff..ab1040c30c 100644
--- a/src/general/semidiscretization.jl
+++ b/src/general/semidiscretization.jl
@@ -690,7 +690,7 @@ end
 function check_update_callback(semi)
     foreach_system(semi) do system
         # This check will be optimized away if the system does not require the callback
-        if requires_update_callback(system) && !semi.update_callback_used[]
+        if requires_update_callback(system, semi) && !semi.update_callback_used[]
             system_name = system |> typeof |> nameof
             throw(ArgumentError("`UpdateCallback` is required for `$system_name`"))
         end
diff --git a/src/preprocessing/particle_packing/system.jl b/src/preprocessing/particle_packing/system.jl
index d6439741ac..65c5160cb2 100644
--- a/src/preprocessing/particle_packing/system.jl
+++ b/src/preprocessing/particle_packing/system.jl
@@ -214,7 +214,7 @@ end
     return ndims(system)
 end
 
-@inline requires_update_callback(system::ParticlePackingSystem) = true
+@inline requires_update_callback(system::ParticlePackingSystem, semi) = true
 
 function write2vtk!(vtk, v, u, t, system::ParticlePackingSystem)
     vtk["velocity"] = [advection_velocity(v, system, particle)
diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl
index d74a69b34b..087f258565 100644
--- a/src/schemes/boundary/open_boundary/system.jl
+++ b/src/schemes/boundary/open_boundary/system.jl
@@ -244,7 +244,7 @@ end
 @inline buffer(system::OpenBoundarySystem) = system.buffer
 
 # The `UpdateCallback` is required to update particle positions between time steps
-@inline requires_update_callback(system::OpenBoundarySystem) = true
+@inline requires_update_callback(system::OpenBoundarySystem, semi) = true
 
 function smoothing_length(system::OpenBoundarySystem, particle)
     return system.smoothing_length
diff --git a/src/schemes/boundary/wall_boundary/dummy_particles.jl b/src/schemes/boundary/wall_boundary/dummy_particles.jl
index 9f5d8926a1..c88658eb28 100644
--- a/src/schemes/boundary/wall_boundary/dummy_particles.jl
+++ b/src/schemes/boundary/wall_boundary/dummy_particles.jl
@@ -681,8 +681,10 @@ end
     (; volume, wall_velocity) = cache
 
     # Prescribed velocity of the boundary particle.
-    # This velocity is zero when not using moving boundaries.
-    v_boundary = current_velocity(v, system, particle)
+    # This velocity is zero when not using moving boundaries or TLSPH.
+    # If not using TLSPH with velocity averaging, this function simply
+    # forwards to `current_velocity`.
+    v_boundary = velocity_for_viscosity(v, system, particle)
 
     for dim in eachindex(v_boundary)
         # The second term is the precalculated smoothed velocity field of the fluid
diff --git a/src/schemes/fluid/shifting_techniques.jl b/src/schemes/fluid/shifting_techniques.jl
index ceaff12a21..3aecb4ea45 100644
--- a/src/schemes/fluid/shifting_techniques.jl
+++ b/src/schemes/fluid/shifting_techniques.jl
@@ -5,7 +5,9 @@ abstract type AbstractShiftingTechnique end
 
 # WARNING: Be careful if defining this function for a specific system type.
 # The version for a specific system type will override this generic version.
-requires_update_callback(system) = requires_update_callback(shifting_technique(system))
+function requires_update_callback(system, semi)
+    return requires_update_callback(shifting_technique(system))
+end
 requires_update_callback(::Nothing) = false
 requires_update_callback(::AbstractShiftingTechnique) = false
 
diff --git a/src/schemes/structure/total_lagrangian_sph/system.jl b/src/schemes/structure/total_lagrangian_sph/system.jl
index 2d2792cd21..05c7c2da7d 100644
--- a/src/schemes/structure/total_lagrangian_sph/system.jl
+++ b/src/schemes/structure/total_lagrangian_sph/system.jl
@@ -7,7 +7,8 @@
                              acceleration=ntuple(_ -> 0.0, NDIMS),
                              penalty_force=nothing, viscosity=nothing,
                              source_terms=nothing, boundary_model=nothing,
-                             self_interaction_nhs=:default)
+                             self_interaction_nhs=:default,
+                             velocity_averaging=nothing)
 
 System for particles of an elastic structure.
 
@@ -61,6 +62,9 @@ See [Total Lagrangian SPH](@ref tlsph) for more details on the method.
                     `transpose_backend` of the [`PrecomputedNeighborhoodSearch`](@ref)
                     is set to `true` when running on GPUs and `false` on CPUs.
                     Alternatively, a user-defined neighborhood search can be passed here.
+- `velocity_averaging`: Velocity averaging technique to be applied on the velocity field
+                    to obtain an averaged velocity to be used in fluid-structure interaction.
+                    See [the docs](@ref velocity_averaging) for details.
 
 !!! note
     If specifying the clamped particles manually (via `n_clamped_particles`),
@@ -83,7 +87,7 @@ See [Total Lagrangian SPH](@ref tlsph) for more details on the method.
     where `beam` and `clamped_particles` are of type [`InitialCondition`](@ref).
 """
 struct TotalLagrangianSPHSystem{BM, NDIMS, ELTYPE <: Real, IC, ARRAY1D, ARRAY2D, ARRAY3D,
-                                YM, PR, LL, LM, K, PF, V, ST, M, IM, NHS,
+                                YM, PR, LL, LM, K, PF, V, ST, M, IM, NHS, VA,
                                 C} <: AbstractStructureSystem{NDIMS}
     initial_condition   :: IC
     initial_coordinates :: ARRAY2D # Array{ELTYPE, 2}: [dimension, particle]
@@ -109,6 +113,7 @@ struct TotalLagrangianSPHSystem{BM, NDIMS, ELTYPE <: Real, IC, ARRAY1D, ARRAY2D,
     clamped_particles_motion :: M
     clamped_particles_moving :: IM
     self_interaction_nhs     :: NHS
+    velocity_averaging       :: VA
     cache                    :: C
 end
 
@@ -121,7 +126,8 @@ function TotalLagrangianSPHSystem(initial_condition, smoothing_kernel, smoothing
                                                       ndims(smoothing_kernel)),
                                   penalty_force=nothing, viscosity=nothing,
                                   source_terms=nothing, boundary_model=nothing,
-                                  self_interaction_nhs=:default)
+                                  self_interaction_nhs=:default,
+                                  velocity_averaging=nothing)
     NDIMS = ndims(initial_condition)
     ELTYPE = eltype(initial_condition)
     n_particles = nparticles(initial_condition)
@@ -186,7 +192,8 @@ function TotalLagrangianSPHSystem(initial_condition, smoothing_kernel, smoothing
     initialize_prescribed_motion!(clamped_particles_motion, initial_condition_sorted,
                                   n_clamped_particles)
 
-    cache = create_cache_tlsph(clamped_particles_motion, initial_condition_sorted)
+    cache = (; create_cache_tlsph(clamped_particles_motion, initial_condition_sorted)...,
+             create_cache_tlsph(velocity_averaging, initial_condition_sorted)...)
 
     return TotalLagrangianSPHSystem(initial_condition_sorted, initial_coordinates,
                                     current_coordinates, mass, correction_matrix,
@@ -197,7 +204,7 @@ function TotalLagrangianSPHSystem(initial_condition, smoothing_kernel, smoothing
                                     smoothing_length, acceleration_, boundary_model,
                                     penalty_force, viscosity, source_terms,
                                     clamped_particles_motion, ismoving,
-                                    self_interaction_nhs, cache)
+                                    self_interaction_nhs, velocity_averaging, cache)
 end
 
 # Initialize self-interaction neighborhood search if not provided by the user
@@ -259,7 +266,8 @@ function initialize_self_interaction_nhs(system::TotalLagrangianSPHSystem,
                                     system.viscosity, system.source_terms,
                                     system.clamped_particles_motion,
                                     system.clamped_particles_moving,
-                                    self_interaction_nhs, system.cache)
+                                    self_interaction_nhs, system.velocity_averaging,
+                                    system.cache)
 end
 
 extract_periodic_box(::Nothing) = nothing
@@ -274,6 +282,22 @@ function create_cache_tlsph(::PrescribedMotion, initial_condition)
     return (; velocity, acceleration)
 end
 
+function create_cache_tlsph(::VelocityAveraging, initial_condition)
+    averaged_velocity = zero(initial_condition.velocity)
+    t_last_averaging = Ref(zero(eltype(initial_condition)))
+
+    return (; averaged_velocity, t_last_averaging)
+end
+
+@inline function requires_update_callback(system::TotalLagrangianSPHSystem, semi)
+    # Velocity averaging used and `integrate_tlsph == false` means that split integration
+    # is used and the averaged velocity is updated there.
+    # Velocity averaging used, `integrate_tlsph == true`, and no fluid system in the semi
+    # means we are inside the split integration
+    return !isnothing(system.velocity_averaging) && semi.integrate_tlsph[] &&
+           any(system -> system isa AbstractFluidSystem, semi.systems)
+end
+
 @inline function Base.eltype(::TotalLagrangianSPHSystem{<:Any, <:Any, ELTYPE}) where {ELTYPE}
     return ELTYPE
 end
@@ -384,6 +408,10 @@ end
     return poisson_ratio[particle]
 end
 
+@inline function velocity_averaging(system::TotalLagrangianSPHSystem)
+    return system.velocity_averaging
+end
+
 function initialize!(system::TotalLagrangianSPHSystem, semi)
     (; correction_matrix) = system
 
diff --git a/src/schemes/structure/total_lagrangian_sph/total_lagrangian_sph.jl b/src/schemes/structure/total_lagrangian_sph/total_lagrangian_sph.jl
index 168d6b634e..a9b859b661 100644
--- a/src/schemes/structure/total_lagrangian_sph/total_lagrangian_sph.jl
+++ b/src/schemes/structure/total_lagrangian_sph/total_lagrangian_sph.jl
@@ -1,5 +1,6 @@
 # Penalty force needs to be included first, so that `dv_penalty_force` is available
 # in the closure of `foreach_point_neighbor`.
 include("penalty_force.jl")
+include("velocity_averaging.jl")
 include("system.jl")
 include("viscosity.jl")
diff --git a/src/schemes/structure/total_lagrangian_sph/velocity_averaging.jl b/src/schemes/structure/total_lagrangian_sph/velocity_averaging.jl
new file mode 100644
index 0000000000..a496ac1ead
--- /dev/null
+++ b/src/schemes/structure/total_lagrangian_sph/velocity_averaging.jl
@@ -0,0 +1,120 @@
+@doc raw"""
+    VelocityAveraging(; time_constant)
+
+Exponential moving average (EMA) of the structure velocity used **only** for the
+fluid-structure viscous coupling (no-slip boundary condition) of a
+[`TotalLagrangianSPHSystem`](@ref).
+
+This is a stabilization option that is typically only needed for challenging FSI
+cases with very stiff structures, where high-frequency noise in the structure velocity
+can trigger instabilities or spurious pressure waves (due to aliasing) in the fluid.
+
+See [the docs](@ref velocity_averaging) for more details.
+
+!!! note "Callback Required"
+    Velocity averaging requires either an [`UpdateCallback`](@ref)
+    or a [`SplitIntegrationCallback`](@ref) to be used in the simulation.
+    In the typical use case of stabilizing a challenging FSI simulation, a
+    [`SplitIntegrationCallback`](@ref) will usually be used anyway for performance reasons.
+
+# Keywords
+- `time_constant`: Time constant ``\tau`` of the EMA (same units as simulation time).
+                   Smaller ``\tau`` smooths only very fast oscillations, which is often
+                   sufficient for stabilizing FSI simulations; larger ``\tau``
+                   can also suppress slower oscillations that are causing spurious
+                   pressure waves without destabilizing the simulation.
+
+# Examples
+```jldoctest; output=false
+VelocityAveraging(time_constant=1e-5)
+
+# output
+VelocityAveraging{Float64}(1.0e-5)
+```
+"""
+struct VelocityAveraging{ELTYPE}
+    time_constant::ELTYPE
+
+    function VelocityAveraging(; time_constant)
+        return new{typeof(time_constant)}(time_constant)
+    end
+end
+
+# No velocity averaging for a system by default
+@inline function velocity_averaging(system)
+    return nothing
+end
+
+@inline function velocity_for_viscosity(v, system, particle)
+    return velocity_for_viscosity(v, velocity_averaging(system), system, particle)
+end
+
+@inline function velocity_for_viscosity(v, ::Nothing, system, particle)
+    return current_velocity(v, system, particle)
+end
+
+@inline function velocity_for_viscosity(v, ::VelocityAveraging, system, particle)
+    if particle <= system.n_integrated_particles
+        return extract_svector(system.cache.averaged_velocity, system, particle)
+    end
+
+    return current_clamped_velocity(v, system, system.clamped_particles_motion, particle)
+end
+
+function initialize_averaged_velocity!(system, v_ode, semi, t)
+    initialize_averaged_velocity!(system, velocity_averaging(system), v_ode, semi, t)
+end
+
+initialize_averaged_velocity!(system, velocity_averaging, v_ode, semi, t) = system
+
+function initialize_averaged_velocity!(system, ::VelocityAveraging, v_ode, semi, t)
+    v = wrap_v(v_ode, system, semi)
+
+    # Make sure the averaged velocity is zero for non-integrated particles
+    set_zero!(system.cache.averaged_velocity)
+
+    @threaded semi for particle in each_integrated_particle(system)
+        v_particle = current_velocity(v, system, particle)
+
+        for i in eachindex(v_particle)
+            system.cache.averaged_velocity[i, particle] = v_particle[i]
+        end
+    end
+    system.cache.t_last_averaging[] = t
+
+    return system
+end
+
+function compute_averaged_velocity!(system, v_ode, semi, t_new)
+    compute_averaged_velocity!(system, velocity_averaging(system), v_ode, semi, t_new)
+end
+
+compute_averaged_velocity!(system, velocity_averaging, v_ode, semi, t_new) = system
+
+function compute_averaged_velocity!(system, velocity_averaging::VelocityAveraging,
+                                    v_ode, semi, t_new)
+    time_constant = velocity_averaging.time_constant
+    averaged_velocity = system.cache.averaged_velocity
+    dt = t_new - system.cache.t_last_averaging[]
+    system.cache.t_last_averaging[] = t_new
+
+    @trixi_timeit timer() "compute averaged velocity" begin
+        v = wrap_v(v_ode, system, semi)
+
+        # Compute an exponential moving average of the velocity with the given time constant
+        alpha = 1 - exp(-dt / time_constant)
+
+        @threaded semi for particle in each_integrated_particle(system)
+            v_particle = @inbounds current_velocity(v, system, particle)
+
+            for i in eachindex(v_particle)
+                @inbounds averaged_velocity[i,
+                                            particle] = (1 - alpha) *
+                                                        averaged_velocity[i, particle] +
+                                                        alpha * v_particle[i]
+            end
+        end
+    end
+
+    return system
+end
diff --git a/test/examples/examples.jl b/test/examples/examples.jl
index c7d9fadb92..55c0707387 100644
--- a/test/examples/examples.jl
+++ b/test/examples/examples.jl
@@ -143,6 +143,8 @@
                                              plate_position=(0.2, 0.0),
                                              tspan=(0.0, 0.2),
                                              E=1e7, # Stiffer plate
+                                             # Test velocity averaging as well
+                                             velocity_averaging=VelocityAveraging(time_constant=1e-4),
                                              maxiters=400,
                                              extra_callback=split_integration) [
                 r"\[ Info: To create the self-interaction neighborhood search.*\n"
@@ -156,6 +158,15 @@
                 @test count_rhs_allocations(sol, semi) == 0
             end
 
+            # Verify that velocity averaging is actually used
+            system = sol.prob.p.systems[3]
+            @test system isa TotalLagrangianSPHSystem
+            # Test that `current_velocity` is not the same as `velocity_for_viscosity`.
+            # `current_velocity` should fail because it tries to access `v`, for which
+            # we passed `nothing`.
+            @test_throws "no method" TrixiParticles.current_velocity(nothing, system, 1)
+            @test !iszero(TrixiParticles.velocity_for_viscosity(nothing, system, 1))
+
             # Now use split integration and verify that it is actually used for TLSPH
             # by using a time step that is too large and verifying that it is crashing.
             split_integration = SplitIntegrationCallback(RDPK3SpFSAL35(), adaptive=false,