Transport Velocity Formalation for WCSPH

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/systems/entropically_damped_sph.md

@@ -54,3 +54,63 @@ Pages = [joinpath("schemes", "fluid", "entropically_damped_sph", "system.jl")]
 - Jonathan R. Clausen. "Entropically damped form of artificial compressibility for explicit simulation of incompressible flow".
   In: American Physical Society 87 (2013), page 13309.
   [doi: 10.1103/PhysRevE.87.013309](http://dx.doi.org/10.1103/PhysRevE.87.013309)
+
+## [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.
+Adami et al (2013) modified the advection velocity and added an extra term in the momentum equation to avoid this problems.
+The authors introduced the so called Transport Velocity Formulation (TVF) for WCSPH. Ramachandran (2019) et al 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 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}
+```
+Note that although ``\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 on the order of the reference pressure but can be chosen arbitrarily large when 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)`` 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]
+```
+where ``V_a``, ``V_b`` denote the volume of particles ``a`` and ``b`` respectively and ``v_{ab}= v_a -v_b``  is the relative velocity.
+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.
+
+The convection tensor is calculated as ``\bm{A} = \rho v (\tilde{v}-v)^T``.
+
+```@autodocs
+Modules = [TrixiParticles]
+Pages = [joinpath("schemes", "fluid", "transport_velocity.jl")]
+```
+
+### References
+- S. Adami, X. Y. Hu, N. A. Adams.
+  "A transport-velocity formulation for smoothed particle hydrodynamics".
+  In: Journal of Computational Physics 241, (2013), pages 292--307.
+  [doi: 10.1016/j.jcp.2013.01.043](http://dx.doi.org/10.1016/j.jcp.2013.01.043)
+- Prabhu Ramachandran. "Entropically damped artificial compressibility for SPH".
+  In: Computers and Fluids 179 (2019), pages 579--594.
+  [doi: 10.1016/j.compfluid.2018.11.023](https://doi.org/10.1016/j.compfluid.2018.11.023)

examples/fluid/lid_driven_cavity_2d.jl

@@ -0,0 +1,113 @@
+# 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
+
+wcsph = true
+
+# ==========================================================================================
+# ==== Resolution
+particle_spacing = 0.02
+
+# 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)
+reynolds_number = 100.0
+
+cavity_size = (1.0, 1.0)
+
+fluid_density = 1.0
+
+const VELOCITY_LID = 1.0
+sound_speed = 10 * VELOCITY_LID
+
+viscosity = ViscosityAdami(; nu=VELOCITY_LID / reynolds_number)
+
+pressure = sound_speed^2 * fluid_density
+
+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}()
+
+if wcsph
+    state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density,
+                                       exponent=1)
+    fluid_system = WeaklyCompressibleSPHSystem(cavity.fluid, ContinuityDensity(),
+                                               state_equation, smoothing_kernel,
+                                               smoothing_length, viscosity=viscosity,
+                                               transport_velocity=TransportVelocityAdami(pressure))
+else
+    state_equation = nothing
+    fluid_system = EntropicallyDampedSPHSystem(cavity.fluid, smoothing_kernel,
+                                               smoothing_length,
+                                               sound_speed, viscosity=viscosity,
+                                               transport_velocity=TransportVelocityAdami(pressure))
+end
+# ==========================================================================================
+# ==== Boundary
+
+movement_function(t) = SVector(VELOCITY_LID * t, 0.0)
+
+is_moving(t) = true
+
+movement = BoundaryMovement(movement_function, is_moving)
+
+boundary_model_cavity = BoundaryModelDummyParticles(cavity.boundary.density,
+                                                    cavity.boundary.mass,
+                                                    AdamiPressureExtrapolation(),
+                                                    viscosity=viscosity,
+                                                    state_equation=state_equation,
+                                                    smoothing_kernel, smoothing_length)
+
+boundary_model_lid = BoundaryModelDummyParticles(lid.density, lid.mass,
+                                                 AdamiPressureExtrapolation(),
+                                                 viscosity=viscosity,
+                                                 state_equation=state_equation,
+                                                 smoothing_kernel, smoothing_length)
+
+boundary_system_cavity = BoundarySPHSystem(cavity.boundary, boundary_model_cavity)
+
+boundary_system_lid = BoundarySPHSystem(lid, boundary_model_lid, movement=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)
+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-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
+            save_everystep=false, callback=callbacks);

examples/fluid/periodic_array_of_cylinders_2d.jl

@@ -0,0 +1,88 @@
+# 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)
+
+const acceleration_x = 2.5e-4
+
+# Boundary geometry and initial fluid particle positions
+cylinder_radius = 0.02
+tank_size = (6cylinder_radius, 4cylinder_radius)
+fluid_size = tank_size
+initial_velocity = (1.2e-4, 0.0)
+
+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);