Add ideal gas equation

examples/fluid/dam_break_2phase_2d.jl

@@ -0,0 +1,89 @@
+# 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_eos = StateEquationCole(; sound_speed, reference_density=air_density, exponent=1,
+                            clip_negative_pressure=false)
+#air_eos = StateEquationIdealGas(; sound_speed, reference_density=air_density, gamma=1.4)
+
+air_system_system = WeaklyCompressibleSPHSystem(air_system, fluid_density_calculator,
+                                                air_eos, 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,15 +17,20 @@ 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)
 
 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, density_diffusion=nothing, sound_speed=sound_speed)
 
 # ==========================================================================================
 # ==== Setup oil layer
@@ -42,8 +47,10 @@ for i in axes(oil.coordinates, 2)
     oil.coordinates[:, i] .+= [0.0, H]
 end
 
+oil_state_equation = StateEquationCole(; sound_speed, reference_density=oil_density,
+                                       exponent=1, clip_negative_pressure=false)
 oil_system = WeaklyCompressibleSPHSystem(oil, fluid_density_calculator,
-                                         state_equation, smoothing_kernel,
+                                         oil_state_equation, smoothing_kernel,
                                          smoothing_length, viscosity=oil_viscosity,
                                          density_diffusion=density_diffusion,
                                          acceleration=(0.0, -gravity),

src/TrixiParticles.jl

@@ -65,7 +65,7 @@ export PenaltyForceGanzenmueller, TransportVelocityAdami
 export SchoenbergCubicSplineKernel, SchoenbergQuarticSplineKernel,
        SchoenbergQuinticSplineKernel, GaussianKernel, WendlandC2Kernel, WendlandC4Kernel,
        WendlandC6Kernel, SpikyKernel, Poly6Kernel
-export StateEquationCole
+export StateEquationCole, StateEquationIdealGas
 export ArtificialViscosityMonaghan, ViscosityAdami, ViscosityMorris
 export DensityDiffusion, DensityDiffusionMolteniColagrossi, DensityDiffusionFerrari,
        DensityDiffusionAntuono

src/schemes/fluid/weakly_compressible_sph/state_equations.jl

@@ -6,9 +6,9 @@ Equation of state to describe the relationship between pressure and density
 of water up to high pressures.
 
 # Keywords
-- `sound_speed`: Artificial speed of sound.
-- `reference_density`: Reference density of the fluid.
-- `exponent`: A value of `7` is usually used for most simulations.
+- `sound_speed`:             Artificial speed of sound.
+- `reference_density`:       Reference density of the fluid.
+- `exponent`:                A value of `7` is usually used for most simulations.
 - `background_pressure=0.0`: Background pressure.
 """
 struct StateEquationCole{ELTYPE, CLIP} # Boolean to clip negative pressure
@@ -49,3 +49,50 @@ function inverse_state_equation(state_equation::StateEquationCole, pressure)
 
     return reference_density * tmp^(1 / exponent)
 end
+
+@doc raw"""
+    StateEquationIdealGas(;sound_speed, reference_density, gamma, background_pressure=0.0,
+                           clip_negative_pressure=false)
+
+Equation of state to describe the relationship between pressure and density
+of a gas using the Ideal Gas Law.
+
+# Keywords
+- `sound_speed`            : Artificial speed of sound.
+- `reference_density`      : Reference density of the fluid.
+- `gamma`                  : Heat-capacity ratio
+- `background_pressure=0.0`: Background pressure.
+"""
+struct StateEquationIdealGas{ELTYPE, CLIP}
+    sound_speed         :: ELTYPE
+    reference_density   :: ELTYPE
+    gamma               :: ELTYPE
+    background_pressure :: ELTYPE
+
+    function StateEquationIdealGas(; sound_speed, reference_density, gamma,
+                                   background_pressure=0.0, clip_negative_pressure=false)
+        new{typeof(sound_speed), clip_negative_pressure}(sound_speed, reference_density,
+                                                         gamma, background_pressure)
+    end
+end
+
+clip_negative_pressure(::StateEquationIdealGas{<:Any, CLIP}) where {CLIP} = CLIP
+
+function (state_equation::StateEquationIdealGas)(density)
+    (; reference_density, sound_speed, gamma, background_pressure) = state_equation
+    pressure = (density - reference_density) * sound_speed^2 / gamma + background_pressure
+
+    # This is determined statically and has therefore no overhead
+    if clip_negative_pressure(state_equation)
+        return max(0.0, pressure)
+    end
+
+    return pressure
+end
+
+function inverse_state_equation(state_equation::StateEquationIdealGas, pressure)
+    (; reference_density, sound_speed, gamma, background_pressure) = state_equation
+    density = (pressure - background_pressure) * gamma / sound_speed^2 + reference_density
+
+    return density
+end

src/visualization/write2vtk.jl

@@ -260,12 +260,19 @@ function write2vtk!(vtk, v, u, t, system::FluidSystem; write_meta_data=true)
             if system.correction isa AkinciFreeSurfaceCorrection
                 vtk["correction_rho0"] = system.correction.rho0
             end
+
+            if system.state_equation isa StateEquationCole
+                vtk["state_equation_exponent"] = system.state_equation.exponent
+            end
+
+            if system.state_equation isa StateEquationIdealGas
+                vtk["state_equation_gamma"] = system.state_equation.gamma
+            end
+
             vtk["state_equation"] = type2string(system.state_equation)
             vtk["state_equation_rho0"] = system.state_equation.reference_density
             vtk["state_equation_pa"] = system.state_equation.background_pressure
             vtk["state_equation_c"] = system.state_equation.sound_speed
-            vtk["state_equation_exponent"] = system.state_equation.exponent
-
             vtk["solver"] = "WCSPH"
         else
             vtk["solver"] = "EDAC"

test/examples/examples.jl

@@ -203,6 +203,18 @@
             @test count_rhs_allocations(sol, semi) == 0
         end
 
+        @trixi_testset "fluid/dam_break_2phase_2d.jl" begin
+            @test_nowarn_mod trixi_include(@__MODULE__,
+                                           joinpath(examples_dir(), "fluid",
+                                                    "dam_break_2phase_2d.jl"),
+                                           tspan=(0.0, 0.05)) [
+                r"┌ Info: The desired tank length in y-direction .*\n",
+                r"└ New tank length in y-direction.*\n"
+            ]
+            @test sol.retcode == ReturnCode.Success
+            @test count_rhs_allocations(sol, semi) == 0
+        end
+
         @trixi_testset "fluid/dam_break_3d.jl" begin
             @test_nowarn_mod trixi_include(@__MODULE__,
                                            joinpath(examples_dir(), "fluid",

test/schemes/fluid/weakly_compressible_sph/state_equation.jl

@@ -139,4 +139,37 @@
             end
         end
     end
+    @testset verbose=true "StateEquationIdealGas" begin
+        # This is the classical ideal gas equation giving a linear relationship
+        # between density and pressure.
+        @testset "Physical Properties of Water" begin
+            # Standard speed of sound for air at 20°C and 1 atm
+            sound_speed = 343.3
+
+            # Density of air at 20°C and 1 atm
+            rest_density = 1.205
+
+            # Heat capacity ratio of air
+            gamma = 1.4
+
+            # Work with pressures in ATM
+            ATM = 101_325.0
+
+            state_equation = StateEquationIdealGas(; sound_speed, gamma=gamma,
+                                                   reference_density=rest_density,
+                                                   background_pressure=1ATM)
+
+            # Results by manual calculation
+            @test TrixiParticles.inverse_state_equation(state_equation, 1ATM) ==
+                  rest_density
+            @test TrixiParticles.inverse_state_equation(state_equation, 100ATM) ==
+                  120.36547777058719
+            @test TrixiParticles.inverse_state_equation(state_equation, 500ATM) ==
+                  601.8219536113436
+
+            @test state_equation(rest_density) == 1ATM
+            @test state_equation(120.36547777058719) == 100ATM
+            @test state_equation(601.8219536113436) == 500ATM
+        end
+    end
 end