Skip to content

Commit

Permalink
Reopen #221: Transport Velocity for EDAC (#436)
Browse files Browse the repository at this point in the history
* add average pressure

* add transport velocity

* add tvf in rhs

* add update callback

* use `step_limiter!` instead of using a callback

* add taylor green vortex example

* add lid driven cavity example

* add periodic array of cylinder example

* add tests

* add Random package

* reexport `seed!`

* remove obsolet `system_index`

* skip CI

* add callback

* rename edac cache

* move functions to `transport_velocity.jl`

* prepare for wcsph

* calculate volume term on the fly

* adapt example files

* time dependent initial velocity function

* time dependent pressure function

* remove velocity function

* multi dimensional functions

* apply formatter

* fix for open boundaries

* add `UpdateCallback`

* fix typo

* prepare for merge `update-callback`

* fix bug

* fix example

* fix tests

* fix update bug

* apply formatter

* fix test

* minor changes

* add tests

* add docs

* fix tpos

* add configuration check

* add setter for tvf

* fix callback used check

* adapt examples

* minor changes

* clean up

* add tgv validation example

* add ldc validation

* modify validation

* increase `maxiters`

* remove double velocity initialization

* add random displacement

* implement suggestions

* update `NEWS.md`

* remove check for viscosity

* fix tests

* add short tilde description in docs

* remove discarded example from tests

* implement suggestions

* apply formatter

* implement suggestions

* add comment

* fix test
  • Loading branch information
LasNikas authored Oct 30, 2024
1 parent 3b4d6fd commit ccd08be
Show file tree
Hide file tree
Showing 23 changed files with 958 additions and 18 deletions.
5 changes: 5 additions & 0 deletions NEWS.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
3 changes: 2 additions & 1 deletion Project.toml
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "TrixiParticles"
uuid = "66699cd8-9c01-4e9d-a059-b96c86d16b3a"
authors = ["erik.faulhaber <44124897+efaulhaber@users.noreply.github.com>"]
version = "0.2.3-dev"
version = "0.2.3"

This comment has been minimized.

Copy link
@svchb

svchb Oct 30, 2024

Collaborator

[deps]
Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
Expand All @@ -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"
Expand Down
14 changes: 14 additions & 0 deletions docs/src/refs.bib
Original file line number Diff line number Diff line change
Expand Up @@ -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},
Expand Down
61 changes: 61 additions & 0 deletions docs/src/systems/entropically_damped_sph.md
Original file line number Diff line number Diff line change
Expand Up @@ -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")]
```
104 changes: 104 additions & 0 deletions examples/fluid/lid_driven_cavity_2d.jl
Original file line number Diff line number Diff line change
@@ -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);
87 changes: 87 additions & 0 deletions examples/fluid/periodic_array_of_cylinders_2d.jl
Original file line number Diff line number Diff line change
@@ -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);
Loading

1 comment on commit ccd08be

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/118371

Tip: Release Notes

Did you know you can add release notes too? Just add markdown formatted text underneath the comment after the text
"Release notes:" and it will be added to the registry PR, and if TagBot is installed it will also be added to the
release that TagBot creates. i.e.

@JuliaRegistrator register

Release notes:

## Breaking changes

- blah

To add them here just re-invoke and the PR will be updated.

Tagging

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v0.2.3 -m "<description of version>" ccd08be91e3e91752ff9afdb8b4280a87820f74d
git push origin v0.2.3

Please sign in to comment.