Add ideal gas equation

examples/fluid/dam_break_2d.jl

@@ -35,7 +35,8 @@ tank_size = (floor(5.366 * H / boundary_particle_spacing) * boundary_particle_sp
 fluid_density = 1000.0
 sound_speed = 20 * sqrt(gravity * H)
 state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density,
-                                   exponent=1, clip_negative_pressure=false)
+                                   exponent=1, clip_negative_pressure=false,
+                                   background_pressure=0.0)
 
 tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fluid_density,
                        n_layers=boundary_layers, spacing_ratio=spacing_ratio,

examples/fluid/dam_break_2phase_2d.jl

@@ -0,0 +1,95 @@
+# 2D dam break simulation with an air layer on top
+
+using TrixiParticles
+using OrdinaryDiffEq
+
+# Size parameters
+H = 0.6
+W = 2 * H
+
+gravity = 9.81
+tspan = (0.0, 2.0)
+
+# Resolution
+# Note: The resolution is very coarse. A better result is obtained with H/60 (which takes over 1 hour)
+fluid_particle_spacing = H / 20
+
+# Numerical settings
+smoothing_length = 3.5 * fluid_particle_spacing
+#sound_speed = 100
+# when using the Ideal gas equation
+sound_speed = 400
+
+# 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, cfl=0.8, exponent=7, background_pressure=0.0,
+              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)
+
+# We use background_pressure here to prevent negative pressure regions in the gas phase.
+air_system_system = WeaklyCompressibleSPHSystem(air_system, fluid_density_calculator,
+                                                StateEquationCole(; sound_speed,
+                                                                  reference_density=air_density,
+                                                                  exponent=1,
+                                                                  clip_negative_pressure=false,
+                                                                  background_pressure=1000),
+                                                smoothing_kernel, smoothing_length,
+                                                viscosity=air_viscosity,
+                                                acceleration=(0.0, -gravity))
+
+# air_system_system = WeaklyCompressibleSPHSystem(air_system, fluid_density_calculator,
+#                                          StateEquationIdealGas(; gamma=1.4,
+#                                                            gas_constant=287.0,
+#                                                            temperature=293.15),
+#                                          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

@@ -64,7 +64,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/density_diffusion.jl

@@ -208,15 +208,15 @@ end
 
     (; delta) = density_diffusion
     (; smoothing_length, state_equation) = particle_system
-    (; sound_speed) = state_equation
+    sound_speed_ = sound_speed(state_equation)
 
     volume_b = m_b / rho_b
 
     psi = density_diffusion_psi(density_diffusion, rho_a, rho_b, pos_diff, distance,
                                 particle_system, particle, neighbor)
     density_diffusion_term = dot(psi, grad_kernel) * volume_b
 
-    dv[end, particle] += delta * smoothing_length * sound_speed * density_diffusion_term
+    dv[end, particle] += delta * smoothing_length * sound_speed_ * density_diffusion_term
 end
 
 # Density diffusion `nothing` or interaction other than fluid-fluid

src/schemes/fluid/weakly_compressible_sph/rhs.jl

@@ -7,7 +7,8 @@ function interact!(dv, v_particle_system, u_particle_system,
                    particle_system::WeaklyCompressibleSPHSystem,
                    neighbor_system)
     (; density_calculator, state_equation, correction, surface_tension) = particle_system
-    (; sound_speed) = state_equation
+
+    sound_speed_ = sound_speed(state_equation)
 
     surface_tension_a = surface_tension_model(particle_system)
     surface_tension_b = surface_tension_model(neighbor_system)
@@ -56,7 +57,7 @@ function interact!(dv, v_particle_system, u_particle_system,
                         dv_viscosity(particle_system, neighbor_system,
                                      v_particle_system, v_neighbor_system,
                                      particle, neighbor, pos_diff, distance,
-                                     sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel)
+                                     sound_speed_, m_a, m_b, rho_a, rho_b, grad_kernel)
 
         dv_surface_tension = surface_tension_correction *
                              surface_tension_force(surface_tension_a, surface_tension_b,

src/schemes/fluid/weakly_compressible_sph/state_equations.jl

@@ -49,3 +49,65 @@ function inverse_state_equation(state_equation::StateEquationCole, pressure)
 
     return reference_density * tmp^(1 / exponent)
 end
+
+@inline sound_speed(eos) = eos.sound_speed
+
+@doc raw"""
+    StateEquationIdealGas(; gas_constant, temperature, gamma)
+Equation of state to describe the relationship between pressure and density
+of a gas using the Ideal Gas Law.
+# Keywords
+- `gas_constant`: Specific gas constant (R) for the gas, typically in J/(kg*K).
+- `temperature` : Absolute temperature of the gas in Kelvin.
+- `gamma`       : Heat-capacity ratio
+This struct calculates the pressure of a gas from its density using the formula:
+\[ P = \rho \cdot R \cdot T \]
+where \( P \) is pressure, \( \rho \) is density, \( R \) is the gas constant, and \( T \) is temperature.
+Note:
+For basic WCSPH this boils down to the assumption of a linear pressure-density relationship.
+"""
+struct StateEquationIdealGas{ELTYPE}
+    gas_constant :: ELTYPE
+    temperature  :: ELTYPE
+    gamma        :: ELTYPE
+
+    function StateEquationIdealGas(; gas_constant, temperature, gamma)
+        new{typeof(gas_constant)}(gas_constant, temperature, gamma)
+    end
+end
+
+function (state_equation::StateEquationIdealGas)(density)
+    (; gas_constant, temperature) = state_equation
+    pressure = density * gas_constant * temperature
+    return pressure
+end
+
+# This version is for simulations that include a temperature.
+function (state_equation::StateEquationIdealGas)(density, internal_energy)
+    (; gamma) = state_equation
+    pressure = (gamma - 1.0) * density * internal_energy
+    return pressure
+end
+
+function inverse_state_equation(state_equation::StateEquationIdealGas, pressure)
+    (; gas_constant, temperature) = state_equation
+    density = pressure / (gas_constant * temperature)
+    return density
+end
+
+# This version is for simulations that include a temperature.
+function inverse_state_equation(state_equation::StateEquationIdealGas, pressure,
+                                internal_energy)
+    gamma = state_equation.gamma
+
+    density = pressure / ((gamma - 1.0) * internal_energy)
+    return density
+end
+
+@inline sound_speed(eos::StateEquationIdealGas) = sqrt(eos.gamma * eos.gas_constant *
+                                                       eos.temperature)
+# This version is for simulations that include a temperature.
+@inline sound_speed(eos::StateEquationIdealGas, pressure, density) = sqrt(eos.gamma *
+                                                                          pressure /
+                                                                          (density *
+                                                                           eos.gas_constant))

src/schemes/fluid/weakly_compressible_sph/system.jl

@@ -215,7 +215,10 @@ end
     return system.pressure[particle]
 end
 
-@inline system_sound_speed(system::WeaklyCompressibleSPHSystem) = system.state_equation.sound_speed
+@inline system_sound_speed(system::WeaklyCompressibleSPHSystem) = sound_speed(system.state_equation)
+@inline system_sound_speed(system::WeaklyCompressibleSPHSystem, pressure, density) = sound_speed(system.state_equation,
+                                                                                                 pressure,
+                                                                                                 density)
 
 function update_quantities!(system::WeaklyCompressibleSPHSystem, v, u,
                             v_ode, u_ode, semi, t)

src/visualization/write2vtk.jl

@@ -261,10 +261,18 @@ function write2vtk!(vtk, v, u, t, system::FluidSystem; write_meta_data=true)
                 vtk["correction_rho0"] = system.correction.rho0
             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
+            if system.state_equation isa StateEquationCole
+                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
+            end
+
+            if system.state_equation isa StateEquationIdealGas
+                vtk["state_equation_gamma"] = system.state_equation.gamma
+                vtk["state_equation_gas_c"] = system.state_equation.gas_constant
+                vtk["state_equation_temperature"] = system.state_equation.temperature
+            end
 
             vtk["solver"] = "WCSPH"
         else

test/examples/examples.jl

@@ -191,6 +191,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",