Reopen #221: Transport Velocity for EDAC

NEWS.md

@@ -3,6 +3,11 @@
 TrixiParticles.jl follows the interpretation of [semantic versioning (semver)](https://julialang.github.io/Pkg.jl/dev/compatibility/#Version-specifier-format-1)
 used in the Julia ecosystem. Notable changes will be documented in this file for human readability.
 
+## Version 0.2.3
+
+### Highlights
+Transport Velocity Formulation (TVF) based on the work of Ramachandran et al. "Entropically damped artificial compressibility for SPH" (2019) was added.
+
 ## Version 0.2.2
 
 ### Highlights

Project.toml

@@ -21,6 +21,7 @@ MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
 PointNeighbors = "1c4d5385-0a27-49de-8e2c-43b175c8985c"
 Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"

docs/src/refs.bib

@@ -13,6 +13,20 @@ @Article{Adami2012
 }
 
 
+@Article{Adami2013,
+  author    = {S. Adami and X.Y. Hu and N.A. Adams},
+  journal   = {Journal of Computational Physics},
+  title     = {A transport-velocity formulation for smoothed particle hydrodynamics},
+  year      = {2013},
+  month     = {may},
+  pages     = {292--307},
+  volume    = {241},
+  doi       = {10.1016/j.jcp.2013.01.043},
+  groups    = {enhancement},
+  publisher = {Elsevier {BV}},
+}
+
+
 @Article{Akinci2012,
   author    = {Akinci, Nadir and Ihmsen, Markus and Akinci, Gizem and Solenthaler, Barbara and Teschner, Matthias},
   journal   = {ACM Transactions on Graphics},

docs/src/systems/entropically_damped_sph.md

@@ -45,3 +45,64 @@ is a good choice for a wide range of Reynolds numbers (0.0125 to 10000).
 Modules = [TrixiParticles]
 Pages = [joinpath("schemes", "fluid", "entropically_damped_sph", "system.jl")]
 ```
+
+## [Transport Velocity Formulation (TVF)](@id transport_velocity_formulation)
+Standard SPH suffers from problems like tensile instability or the creation of void regions in the flow.
+To address these problems, [Adami (2013)](@cite Adami2013) modified the advection velocity and added an extra term to the momentum equation.
+The authors introduced the so-called Transport Velocity Formulation (TVF) for WCSPH. [Ramachandran (2019)](@cite Ramachandran2019) applied the TVF
+also for the [EDAC](@ref edac) scheme.
+
+The transport velocity ``\tilde{v}_a`` of particle ``a`` is used to evolve the position of the particle ``r_a`` from one time step to the next by
+
+```math
+\frac{\mathrm{d} r_a}{\mathrm{d}t} = \tilde{v}_a
+```
+
+and is obtained at every time-step ``\Delta t`` from
+
+```math
+\tilde{v}_a (t + \Delta t) = v_a (t) + \Delta t \left(\frac{\tilde{\mathrm{d}} v_a}{\mathrm{d}t} - \frac{1}{\rho_a} \nabla p_{\text{background}} \right),
+```
+
+where ``\rho_a`` is the density of particle ``a`` and ``p_{\text{background}}`` is a constant background pressure field.
+The tilde in the second term of the right hand side indicates that the material derivative has an advection part.
+
+The discretized form of the last term is
+
+```math
+ -\frac{1}{\rho_a} \nabla p_{\text{background}} \approx  -\frac{p_{\text{background}}}{m_a} \sum_b \left(V_a^2 + V_b^2 \right) \nabla_a W_{ab},
+```
+
+where ``V_a``, ``V_b`` denote the volume of particles ``a`` and ``b`` respectively.
+Note that although in the continuous case ``\nabla p_{\text{background}} = 0``, the discretization is not 0th-order consistent for **non**-uniform particle distribution,
+which means that there is a non-vanishing contribution only when particles are disordered.
+That also means that ``p_{\text{background}}`` occurs as prefactor to correct the trajectory of a particle resulting in uniform pressure distributions.
+Suggested is a background pressure which is in the order of the reference pressure but can be chosen arbitrarily large when the time-step criterion is adjusted.
+
+The inviscid momentum equation with an additional convection term for a particle moving with ``\tilde{v}`` is
+
+```math
+\frac{\tilde{\mathrm{d}} \left( \rho v \right)}{\mathrm{d}t} = -\nabla p +  \nabla \cdot \bm{A},
+```
+
+ where the tensor ``\bm{A} = \rho v\left(\tilde{v}-v\right)^T`` is a consequence of the modified
+ advection velocity and can be interpreted as the convection of momentum with the relative velocity ``\tilde{v}-v``.
+
+The discretized form of the momentum equation for a particle ``a`` reads as
+
+```math
+\frac{\tilde{\mathrm{d}} v_a}{\mathrm{d}t} = \frac{1}{m_a} \sum_b \left(V_a^2 + V_b^2 \right) \left[ -\tilde{p}_{ab} \nabla_a W_{ab} + \frac{1}{2} \left(\bm{A}_a + \bm{A}_b \right) \cdot \nabla_a W_{ab} \right].
+```
+
+Here, ``\tilde{p}_{ab}`` is the density-weighted pressure
+
+```math
+\tilde{p}_{ab} = \frac{\rho_b p_a + \rho_a p_b}{\rho_a + \rho_b},
+```
+
+with the density  ``\rho_a``,  ``\rho_b`` and the pressure  ``p_a``,  ``p_b`` of particles ``a`` and ``b`` respectively. ``\bm{A}_a`` and ``\bm{A}_b`` are the convection tensors for particle ``a`` and ``b`` respectively and are given, e.g. for particle ``a``, as ``\bm{A}_a = \rho v_a\left(\tilde{v}_a-v_a\right)^T``.
+
+```@autodocs
+Modules = [TrixiParticles]
+Pages = [joinpath("schemes", "fluid", "transport_velocity.jl")]
+```

examples/fluid/lid_driven_cavity_2d.jl

@@ -0,0 +1,104 @@
+# Lid-driven cavity
+#
+# S. Adami et al
+# "A transport-velocity formulation for smoothed particle hydrodynamics".
+# In: Journal of Computational Physics, Volume 241 (2013), pages 292-307.
+# https://doi.org/10.1016/j.jcp.2013.01.043
+
+using TrixiParticles
+using OrdinaryDiffEq
+
+# ==========================================================================================
+# ==== Resolution
+particle_spacing = 0.02
+
+# Make sure that the kernel support of fluid particles at a boundary is always fully sampled
+boundary_layers = 4
+
+# ==========================================================================================
+# ==== Experiment Setup
+tspan = (0.0, 5.0)
+reynolds_number = 100.0
+
+cavity_size = (1.0, 1.0)
+
+fluid_density = 1.0
+
+const VELOCITY_LID = 1.0
+sound_speed = 10 * VELOCITY_LID
+
+pressure = sound_speed^2 * fluid_density
+
+viscosity = ViscosityAdami(; nu=VELOCITY_LID / reynolds_number)
+
+cavity = RectangularTank(particle_spacing, cavity_size, cavity_size, fluid_density,
+                         n_layers=boundary_layers,
+                         faces=(true, true, true, false), pressure=pressure)
+
+lid_position = 0.0 - particle_spacing * boundary_layers
+lid_length = cavity.n_particles_per_dimension[1] + 2boundary_layers
+
+lid = RectangularShape(particle_spacing, (lid_length, 3),
+                       (lid_position, cavity_size[2]), density=fluid_density)
+
+# ==========================================================================================
+# ==== Fluid
+
+smoothing_length = 1.0 * particle_spacing
+smoothing_kernel = SchoenbergQuinticSplineKernel{2}()
+
+fluid_system = EntropicallyDampedSPHSystem(cavity.fluid, smoothing_kernel, smoothing_length,
+                                           density_calculator=ContinuityDensity(),
+                                           sound_speed, viscosity=viscosity,
+                                           transport_velocity=TransportVelocityAdami(pressure))
+
+# ==========================================================================================
+# ==== Boundary
+
+lid_movement_function(t) = SVector(VELOCITY_LID * t, 0.0)
+
+is_moving(t) = true
+
+lid_movement = BoundaryMovement(lid_movement_function, is_moving)
+
+boundary_model_cavity = BoundaryModelDummyParticles(cavity.boundary.density,
+                                                    cavity.boundary.mass,
+                                                    AdamiPressureExtrapolation(),
+                                                    viscosity=viscosity,
+                                                    smoothing_kernel, smoothing_length)
+
+boundary_model_lid = BoundaryModelDummyParticles(lid.density, lid.mass,
+                                                 AdamiPressureExtrapolation(),
+                                                 viscosity=viscosity,
+                                                 smoothing_kernel, smoothing_length)
+
+boundary_system_cavity = BoundarySPHSystem(cavity.boundary, boundary_model_cavity)
+
+boundary_system_lid = BoundarySPHSystem(lid, boundary_model_lid, movement=lid_movement)
+
+# ==========================================================================================
+# ==== Simulation
+bnd_thickness = boundary_layers * particle_spacing
+periodic_box = PeriodicBox(min_corner=[-bnd_thickness, -bnd_thickness],
+                           max_corner=cavity_size .+ [bnd_thickness, bnd_thickness])
+
+semi = Semidiscretization(fluid_system, boundary_system_cavity, boundary_system_lid,
+                          neighborhood_search=GridNeighborhoodSearch{2}(; periodic_box))
+
+ode = semidiscretize(semi, tspan)
+
+info_callback = InfoCallback(interval=100)
+
+saving_callback = SolutionSavingCallback(dt=0.02)
+
+pp_callback = nothing
+
+callbacks = CallbackSet(info_callback, saving_callback, pp_callback, UpdateCallback())
+
+# Use a Runge-Kutta method with automatic (error based) time step size control
+sol = solve(ode, RDPK3SpFSAL49(),
+            abstol=1e-6, # Default abstol is 1e-6 (may needs to be tuned to prevent boundary penetration)
+            reltol=1e-4, # Default reltol is 1e-3 (may needs to be tuned to prevent boundary penetration)
+            dtmax=1e-2, # Limit stepsize to prevent crashing
+            maxiters=Int(1e7),
+            save_everystep=false, callback=callbacks);

examples/fluid/periodic_array_of_cylinders_2d.jl

@@ -0,0 +1,87 @@
+# Channel flow through periodic array of cylinders
+#
+# S. Adami et al.
+# "A transport-velocity formulation for smoothed particle hydrodynamics".
+# In: Journal of Computational Physics, Volume 241 (2013), pages 292-307.
+# https://doi.org/10.1016/j.jcp.2013.01.043
+
+using TrixiParticles
+using OrdinaryDiffEq
+
+# ==========================================================================================
+# ==== Resolution
+n_particles_x = 144
+
+# Make sure that the kernel support of fluid particles at a boundary is always fully sampled
+boundary_layers = 3
+
+# ==========================================================================================
+# ==== Experiment Setup
+tspan = (0.0, 5.0)
+
+acceleration_x = 2.5e-4
+
+# Boundary geometry and initial fluid particle positions
+cylinder_radius = 0.02
+tank_size = (6 * cylinder_radius, 4 * cylinder_radius)
+fluid_size = tank_size
+
+fluid_density = 1000.0
+nu = 0.1 / fluid_density # viscosity parameter
+
+# Adami uses `c = 0.1 * sqrt(acceleration_x * cylinder_radius)``  but the original setup
+# from M. Ellero and N. A. Adams (https://doi.org/10.1002/nme.3088) uses `c = 0.02`
+sound_speed = 0.02
+
+pressure = sound_speed^2 * fluid_density
+
+particle_spacing = tank_size[1] / n_particles_x
+
+box = RectangularTank(particle_spacing, fluid_size, tank_size,
+                      fluid_density, n_layers=boundary_layers,
+                      pressure=pressure, faces=(false, false, true, true))
+
+cylinder = SphereShape(particle_spacing, cylinder_radius, tank_size ./ 2,
+                       fluid_density, sphere_type=RoundSphere())
+
+fluid = setdiff(box.fluid, cylinder)
+boundary = union(cylinder, box.boundary)
+
+# ==========================================================================================
+# ==== Fluid
+smoothing_length = 1.2 * particle_spacing
+smoothing_kernel = SchoenbergQuarticSplineKernel{2}()
+fluid_system = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, smoothing_length,
+                                           sound_speed, viscosity=ViscosityAdami(; nu),
+                                           transport_velocity=TransportVelocityAdami(pressure),
+                                           acceleration=(acceleration_x, 0.0))
+
+# ==========================================================================================
+# ==== Boundary
+boundary_model = BoundaryModelDummyParticles(boundary.density, boundary.mass,
+                                             AdamiPressureExtrapolation(),
+                                             viscosity=ViscosityAdami(; nu),
+                                             smoothing_kernel, smoothing_length)
+
+boundary_system = BoundarySPHSystem(boundary, boundary_model)
+
+# ==========================================================================================
+# ==== Simulation
+periodic_box = PeriodicBox(min_corner=[0.0, -tank_size[2]],
+                           max_corner=[tank_size[1], 2 * tank_size[2]])
+semi = Semidiscretization(fluid_system, boundary_system,
+                          neighborhood_search=GridNeighborhoodSearch{2}(; periodic_box))
+
+ode = semidiscretize(semi, tspan)
+
+info_callback = InfoCallback(interval=10)
+saving_callback = SolutionSavingCallback(dt=0.02, prefix="")
+
+callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback())
+
+# Use a Runge-Kutta method with automatic (error based) time step size control
+sol = solve(ode, RDPK3SpFSAL49(),
+            abstol=1e-8, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration)
+            reltol=1e-4, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration)
+            dtmax=1e-2, # Limit stepsize to prevent crashing
+            save_everystep=false, callback=callbacks);

examples/fluid/taylor_green_vortex_2d.jl

@@ -0,0 +1,93 @@
+# Taylor Green vortex
+#
+# P. Ramachandran, K. Puri
+# "Entropically damped artificial compressibility for SPH".
+# In: Computers and Fluids, Volume 179 (2019), pages 579-594.
+# https://doi.org/10.1016/j.compfluid.2018.11.023
+
+using TrixiParticles
+using OrdinaryDiffEq
+
+# ==========================================================================================
+# ==== Resolution
+particle_spacing = 0.02
+
+# ==========================================================================================
+# ==== Experiment Setup
+tspan = (0.0, 5.0)
+reynolds_number = 100.0
+
+box_length = 1.0
+
+U = 1.0 # m/s
+fluid_density = 1.0
+sound_speed = 10U
+
+b = -8pi^2 / reynolds_number
+
+# Taylor Green Vortex Pressure Function
+function pressure_function(pos, t)
+    x = pos[1]
+    y = pos[2]
+
+    return -U^2 * exp(2 * b * t) * (cos(4pi * x) + cos(4pi * y)) / 4
+end
+
+initial_pressure_function(pos) = pressure_function(pos, 0.0)
+
+# Taylor Green Vortex Velocity Function
+function velocity_function(pos, t)
+    x = pos[1]
+    y = pos[2]
+
+    vel = U * exp(b * t) * [-cos(2pi * x) * sin(2pi * y), sin(2pi * x) * cos(2pi * y)]
+
+    return SVector{2}(vel)
+end
+
+initial_velocity_function(pos) = velocity_function(pos, 0.0)
+
+n_particles_xy = round(Int, box_length / particle_spacing)
+
+# ==========================================================================================
+# ==== Fluid
+nu = U * box_length / reynolds_number
+
+background_pressure = sound_speed^2 * fluid_density
+
+smoothing_length = 1.0 * particle_spacing
+smoothing_kernel = SchoenbergQuinticSplineKernel{2}()
+
+fluid = RectangularShape(particle_spacing, (n_particles_xy, n_particles_xy), (0.0, 0.0),
+                         coordinates_perturbation=0.2, # To avoid stagnant streamlines when not using TVF.
+                         density=fluid_density, pressure=initial_pressure_function,
+                         velocity=initial_velocity_function)
+
+fluid_system = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, smoothing_length,
+                                           sound_speed,
+                                           transport_velocity=TransportVelocityAdami(background_pressure),
+                                           viscosity=ViscosityAdami(; nu))
+
+# ==========================================================================================
+# ==== Simulation
+periodic_box = PeriodicBox(min_corner=[0.0, 0.0], max_corner=[box_length, box_length])
+semi = Semidiscretization(fluid_system,
+                          neighborhood_search=GridNeighborhoodSearch{2}(; periodic_box))
+
+ode = semidiscretize(semi, tspan)
+
+info_callback = InfoCallback(interval=100)
+
+saving_callback = SolutionSavingCallback(dt=0.02)
+
+pp_callback = nothing
+
+callbacks = CallbackSet(info_callback, saving_callback, pp_callback, UpdateCallback())
+
+dt_max = min(smoothing_length / 4 * (sound_speed + U), smoothing_length^2 / (8 * nu))
+
+# Use a Runge-Kutta method with automatic (error based) time step size control
+sol = solve(ode, RDPK3SpFSAL49(),
+            abstol=1e-8, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration)
+            reltol=1e-4, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration)
+            dtmax=dt_max, save_everystep=false, callback=callbacks);