Implement Morris surface tension model

NEWS.md

@@ -3,51 +3,83 @@
 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.4
+
+### Features
+
+- Support for surface tension was added to EDAC (#539)
+- A method to prevent penetration of fast moving particles with solids was added (#498) 
+
+### Documentation
+
+- GPU Support Documentation was added (#660)
+
+### Fixes
+
+- JuliaFormatter was updated from 1.0.42 to 1.0.62 (#661)
+- Diverse Doc fixes (#663, #659, #637, #658, #664)
+- Simulations can be run with `Float32` (#662)
+
+### Refactored
+
+- Surface normal calculation was moved from surface_tension.jl to surface_normal_sph.jl (#539)
+
 ## 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
+
 Hotfix for threaded sampling of complex geometries.
 
 ## Version 0.2.1
 
 ### Highlights
+
 Particle sampling of complex geometries from `.stl` and `.asc` files.
 
 ## Version 0.2.0
 
 ### Removed
+
 Use of the internal neighborhood search has been removed and replaced with PointNeighbors.jl.
 
 ## Development Cycle 0.1
 
 ### Highlights
 
 #### Discrete Element Method
+
 A basic implementation of the discrete element method was added.
 
 #### Surface Tension and Adhesion Model
+
 A surface tension and adhesion model based on the work by Akinci et al., "Versatile Surface Tension and Adhesion for SPH Fluids" (2013) was added to WCSPH.
 
 #### Support for Open Boundaries
+
 Open boundaries using the method of characteristics based on the work of Lastiwka et al., "Permeable and non-reflecting boundary conditions in SPH" (2009) were added for WCSPH and EDAC.
 
 ## Pre Initial Release (v0.1.0)
+
 This section summarizes the initial features that TrixiParticles.jl was released with.
 
 ### Highlights
 #### EDAC
+
 An implementation of EDAC (Entropically Damped Artificial Compressibility) was added,
 which allows for more stable simulations compared to basic WCSPH and reduces spurious pressure oscillations.
 
 #### WCSPH
+
 An implementation of WCSPH (Weakly Compressible Smoothed Particle Hydrodynamics), which is the classical SPH approach.
 
 Features:
+
 - Correction schemes (Shepard (0. Order) ... MixedKernelGradient (1. Order))
 - Density reinitialization
 - Kernel summation and Continuity equation density formulations
@@ -57,4 +89,5 @@ Features:
 
 
 #### TLSPH
+
 An implementation of TLSPH (Total Lagrangian Smoothed Particle Hydrodynamics) for solid bodies enabling FSI (Fluid Structure Interactions).

Project.toml

@@ -1,7 +1,7 @@
 name = "TrixiParticles"
 uuid = "66699cd8-9c01-4e9d-a059-b96c86d16b3a"
 authors = ["erik.faulhaber <[email protected]>"]
-version = "0.2.4-dev"
+version = "0.2.4"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"

examples/fluid/dam_break_2d.jl

@@ -61,7 +61,8 @@ fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator,
                                            smoothing_length, viscosity=viscosity,
                                            density_diffusion=density_diffusion,
                                            acceleration=(0.0, -gravity), correction=nothing,
-                                           surface_tension=nothing)
+                                           surface_tension=nothing,
+                                           reference_particle_spacing=fluid_particle_spacing)
 
 # ==========================================================================================
 # ==== Boundary

examples/fluid/dam_break_2phase_2d.jl

@@ -0,0 +1,96 @@
+# 2D dam break simulation with an air layer on top
+
+using TrixiParticles
+using OrdinaryDiffEq
+
+# Size parameters
+H = 0.6
+W = 2 * H
+
+# ==========================================================================================
+# ==== Resolution
+
+# Note: The resolution is very coarse. A better result is obtained with H/60 or higher (which takes over 1 hour+)
+fluid_particle_spacing = H / 20
+
+# ==========================================================================================
+# ==== Experiment Setup
+gravity = 9.81
+tspan = (0.0, 2.0)
+
+# Numerical settings
+smoothing_length = 3.5 * fluid_particle_spacing
+sound_speed = 100.0
+# when using the Ideal gas equation
+# sound_speed = 343.0
+
+# physical values
+nu_water = 8.9E-7
+nu_air = 1.544E-5
+nu_ratio = nu_water / nu_air
+
+nu_sim_air = 0.02 * smoothing_length * sound_speed
+nu_sim_water = nu_ratio * nu_sim_air
+
+air_viscosity = ViscosityMorris(nu=nu_sim_air)
+water_viscosity = ViscosityMorris(nu=nu_sim_water)
+
+trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "dam_break_2d.jl"),
+              sol=nothing, fluid_particle_spacing=fluid_particle_spacing,
+              viscosity=water_viscosity, smoothing_length=smoothing_length,
+              gravity=gravity, tspan=tspan, density_diffusion=nothing,
+              sound_speed=sound_speed, exponent=7,
+              tank_size=(floor(5.366 * H / fluid_particle_spacing) * fluid_particle_spacing,
+                         2.6 * H))
+
+# ==========================================================================================
+# ==== Setup air_system layer
+
+air_size = (tank_size[1], 1.5 * H)
+air_size2 = (tank_size[1] - W, H)
+air_density = 1.0
+
+# Air above the initial water
+air_system = RectangularShape(fluid_particle_spacing,
+                              round.(Int, air_size ./ fluid_particle_spacing),
+                              zeros(length(air_size)), density=air_density)
+
+# Air for the rest of the empty volume
+air_system2 = RectangularShape(fluid_particle_spacing,
+                               round.(Int, air_size2 ./ fluid_particle_spacing),
+                               (W, 0.0), density=air_density)
+
+# move on top of the water
+for i in axes(air_system.coordinates, 2)
+    air_system.coordinates[:, i] .+= [0.0, H]
+end
+
+air_system = union(air_system, air_system2)
+
+air_system_system = WeaklyCompressibleSPHSystem(air_system, fluid_density_calculator,
+                                                StateEquationCole(; sound_speed,
+                                                                  reference_density=air_density,
+                                                                  exponent=1,
+                                                                  clip_negative_pressure=false,
+                                                                  background_pressure=0),
+                                                smoothing_kernel, smoothing_length,
+                                                viscosity=air_viscosity,
+                                                acceleration=(0.0, -gravity))
+
+# air_system_system = WeaklyCompressibleSPHSystem(air_system, fluid_density_calculator,
+#                                          StateEquationIdealGas(; sound_speed, reference_density=air_density, gamma=1.4, background_pressure=0),
+#                                          smoothing_kernel, smoothing_length,
+#                                          viscosity=air_viscosity,
+#                                          acceleration=(0.0, -gravity))
+
+# ==========================================================================================
+# ==== Simulation
+semi = Semidiscretization(fluid_system, air_system_system, boundary_system,
+                          neighborhood_search=GridNeighborhoodSearch{2}(update_strategy=nothing))
+ode = semidiscretize(semi, tspan, data_type=nothing)
+
+sol = solve(ode, RDPK3SpFSAL35(),
+            abstol=1e-5, # 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/dam_break_oil_film_2d.jl

@@ -17,21 +17,39 @@ fluid_particle_spacing = H / 60
 smoothing_length = 3.5 * fluid_particle_spacing
 sound_speed = 20 * sqrt(gravity * H)
 
-nu = 0.02 * smoothing_length * sound_speed / 8
-oil_viscosity = ViscosityMorris(nu=10 * nu)
+# Physical values
+nu_water = 8.9E-7
+nu_oil = 6E-5
+nu_ratio = nu_water / nu_oil
 
+nu_sim_oil = max(0.01 * smoothing_length * sound_speed, nu_oil)
+nu_sim_water = nu_ratio * nu_sim_oil
+
+oil_viscosity = ViscosityMorris(nu=nu_sim_oil)
+# Physical values
+nu_water = 8.9E-7
+nu_oil = 6E-5
+nu_ratio = nu_water / nu_oil
+
+nu_sim_oil = max(0.01 * smoothing_length * sound_speed, nu_oil)
+nu_sim_water = nu_ratio * nu_sim_oil
+
+oil_viscosity = ViscosityMorris(nu=nu_sim_oil)
+
+# TODO: broken if both are set to surface tension
 trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "dam_break_2d.jl"),
               sol=nothing, fluid_particle_spacing=fluid_particle_spacing,
-              viscosity=ViscosityMorris(nu=nu), smoothing_length=smoothing_length,
-              gravity=gravity,
-              density_diffusion=DensityDiffusionMolteniColagrossi(delta=0.1),
-              sound_speed=sound_speed)
+              viscosity=ViscosityMorris(nu=nu_sim_water), smoothing_length=smoothing_length,
+              gravity=gravity, tspan=tspan, density_diffusion=nothing,
+              sound_speed=sound_speed, prefix="")
 
 # ==========================================================================================
 # ==== Setup oil layer
 
 oil_size = (W, 0.1 * H)
 oil_density = 700.0
+oil_eos = StateEquationCole(; sound_speed, reference_density=oil_density, exponent=1,
+                            clip_negative_pressure=false)
 
 oil = RectangularShape(fluid_particle_spacing,
                        round.(Int, oil_size ./ fluid_particle_spacing),
@@ -43,12 +61,19 @@ for i in axes(oil.coordinates, 2)
 end
 
 oil_system = WeaklyCompressibleSPHSystem(oil, fluid_density_calculator,
-                                         state_equation, smoothing_kernel,
+                                         oil_eos, smoothing_kernel,
                                          smoothing_length, viscosity=oil_viscosity,
-                                         density_diffusion=density_diffusion,
                                          acceleration=(0.0, -gravity),
-                                         surface_tension=SurfaceTensionAkinci(surface_tension_coefficient=0.02),
-                                         correction=AkinciFreeSurfaceCorrection(oil_density))
+                                         surface_tension=SurfaceTensionAkinci(surface_tension_coefficient=0.01),
+                                         correction=AkinciFreeSurfaceCorrection(oil_density),
+                                         reference_particle_spacing=fluid_particle_spacing)
+
+# oil_system = WeaklyCompressibleSPHSystem(oil, fluid_density_calculator,
+#                                          oil_eos, smoothing_kernel,
+#                                          smoothing_length, viscosity=oil_viscosity,
+#                                          acceleration=(0.0, -gravity),
+#                                          surface_tension=SurfaceTensionMorris(surface_tension_coefficient=0.03),
+#                                          reference_particle_spacing=fluid_particle_spacing)
 
 # ==========================================================================================
 # ==== Simulation

examples/fluid/falling_water_spheres_2d.jl

@@ -40,7 +40,7 @@ sphere2 = SphereShape(fluid_particle_spacing, sphere_radius, sphere2_center,
 
 # ==========================================================================================
 # ==== Fluid
-fluid_smoothing_length = 1.0 * fluid_particle_spacing
+fluid_smoothing_length = 1.0 * fluid_particle_spacing - eps()
 fluid_smoothing_kernel = SchoenbergCubicSplineKernel{2}()
 
 fluid_density_calculator = ContinuityDensity()
@@ -50,13 +50,13 @@ alpha = 8 * nu / (fluid_smoothing_length * sound_speed)
 viscosity = ArtificialViscosityMonaghan(alpha=alpha, beta=0.0)
 density_diffusion = DensityDiffusionAntuono(sphere2, delta=0.1)
 
-sphere_surface_tension = WeaklyCompressibleSPHSystem(sphere1, fluid_density_calculator,
-                                                     state_equation, fluid_smoothing_kernel,
+sphere_surface_tension = EntropicallyDampedSPHSystem(sphere1, fluid_smoothing_kernel,
                                                      fluid_smoothing_length,
-                                                     viscosity=viscosity,
+                                                     sound_speed, viscosity=viscosity,
+                                                     density_calculator=ContinuityDensity(),
                                                      acceleration=(0.0, -gravity),
                                                      surface_tension=SurfaceTensionAkinci(surface_tension_coefficient=0.05),
-                                                     correction=AkinciFreeSurfaceCorrection(fluid_density))
+                                                     reference_particle_spacing=fluid_particle_spacing)
 
 sphere = WeaklyCompressibleSPHSystem(sphere2, fluid_density_calculator,
                                      state_equation, fluid_smoothing_kernel,
@@ -75,16 +75,16 @@ boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundar
                                              viscosity=ViscosityAdami(nu=wall_viscosity))
 
 boundary_system = BoundarySPHSystem(tank.boundary, boundary_model,
-                                    adhesion_coefficient=0.001)
+                                    adhesion_coefficient=1.0)
 
 # ==========================================================================================
 # ==== Simulation
-semi = Semidiscretization(boundary_system, sphere_surface_tension, sphere)
+semi = Semidiscretization(sphere_surface_tension, sphere, boundary_system)
 ode = semidiscretize(semi, tspan)
 
 info_callback = InfoCallback(interval=50)
-saving_callback = SolutionSavingCallback(dt=0.01, output_directory="out", prefix="",
-                                         write_meta_data=true)
+saving_callback = SolutionSavingCallback(dt=0.01, output_directory="out",
+                                         prefix="", write_meta_data=true)
 
 callbacks = CallbackSet(info_callback, saving_callback)
 

examples/fluid/falling_water_spheres_3d.jl

@@ -3,35 +3,35 @@ using OrdinaryDiffEq
 
 # ==========================================================================================
 # ==== Resolution
-fluid_particle_spacing = 0.008
+fluid_particle_spacing = 0.005
 
 # ==========================================================================================
 # ==== Experiment Setup
 gravity = 9.81
-nu = 0.01
+nu = 0.001
 fluid_density = 1000.0
 sound_speed = 50
 
 sphere1_radius = 0.05
 
-sphere1_center = (0.5, 0.5, 0.2)
-sphere2_center = (1.5, 0.5, 0.2)
+sphere1_center = (0.5, 0.5, 0.075)
+sphere2_center = (1.5, 0.5, 0.075)
 sphere1 = SphereShape(fluid_particle_spacing, sphere1_radius, sphere1_center,
-                      fluid_density, sphere_type=VoxelSphere(), velocity=(0.0, 0.0, -2.0))
+                      fluid_density, sphere_type=VoxelSphere(), velocity=(0.0, 0.0, -1.0))
 sphere2 = SphereShape(fluid_particle_spacing, sphere1_radius, sphere2_center,
-                      fluid_density, sphere_type=VoxelSphere(), velocity=(0.0, 0.0, -2.0))
+                      fluid_density, sphere_type=VoxelSphere(), velocity=(0.0, 0.0, -1.0))
 
 # `compact_support` needs to be `2.0 * particle_spacing` to be correct
 fluid_smoothing_length = 1.0 * fluid_particle_spacing
 
 trixi_include(@__MODULE__,
               joinpath(examples_dir(), "fluid", "falling_water_spheres_2d.jl"),
-              fluid_particle_spacing=fluid_particle_spacing, tspan=(0.0, 0.2),
+              fluid_particle_spacing=fluid_particle_spacing, tspan=(0.0, 0.1),
               initial_fluid_size=(0.0, 0.0, 0.0),
               tank_size=(2.0, 1.0, 0.1), sound_speed=sound_speed,
               faces=(true, true, true, true, true, false),
               acceleration=(0.0, 0.0, -gravity), sphere1=sphere1, sphere2=sphere2,
               fluid_smoothing_length=fluid_smoothing_length,
               fluid_smoothing_kernel=SchoenbergCubicSplineKernel{3}(),
               nu=nu, alpha=10 * nu / (fluid_smoothing_length * sound_speed),
-              surface_tension_coefficient=1.5, adhesion_coefficient=0.5)
+              surface_tension_coefficient=10, adhesion_coefficient=0.1)