diff --git a/examples/fluid/vortex_street_2d.jl b/examples/fluid/vortex_street_2d.jl new file mode 100644 index 0000000000..93ff39a752 --- /dev/null +++ b/examples/fluid/vortex_street_2d.jl @@ -0,0 +1,164 @@ +# ========================================================================================== +# 2D Vortex Street +# +# Based on: +# A. Tafuni, J. M. Domínguez, R. Vacondio, A. J. C. Crespo. +# "A versatile algorithm for the treatment of open boundary conditions in smoothed +# particle hydrodynamics GPU models". +# Computer Methods in Applied Mechanics and Engineering, Volume 342 (2018), pages 604-624. +# https://doi.org/10.1016/j.cma.2018.08.004 +# ========================================================================================== + +using TrixiParticles +using OrdinaryDiffEqLowStorageRK + +# ========================================================================================== +# ==== Resolution +particle_spacing_factor = 0.1 # Resolution in the paper is `0.01` (5M particles) + +cylinder_diameter = 0.1 +particle_spacing = particle_spacing_factor * cylinder_diameter + +# Make sure that the kernel support of fluid particles at a boundary is always fully sampled +boundary_layers = 4 + +# Make sure that the kernel support of fluid particles at an open boundary is always +# fully sampled. +# Note: Due to the dynamics at the inlets and outlets of open boundaries, +# it is recommended to use `open_boundary_layers > boundary_layers`. +open_boundary_layers = 8 + +# ========================================================================================== +# ==== Experiment Setup +tspan = (0.0, 4.0) + +# We use a smaller domain for a faster runtime. The original domain size in the paper is +# domain_size = (25 * cylinder_diameter, 20 * cylinder_diameter) +domain_size = (20 * cylinder_diameter, 10 * cylinder_diameter) +open_boundary_size = (particle_spacing * open_boundary_layers, domain_size[2]) + +flow_direction = SVector(1.0, 0.0) +reynolds_number = 200 +prescribed_velocity = 1.0 +fluid_density = 1000.0 +# Maximum velocity observed in the vortex street is typically around 1.5 +v_max = 1.5 +sound_speed = 10 * v_max + +pipe = RectangularTank(particle_spacing, domain_size, domain_size, fluid_density; + n_layers=boundary_layers, + velocity=prescribed_velocity * flow_direction, + faces=(false, false, true, true), + coordinates_eltype=Float64) + +min_coords_inlet = (-open_boundary_layers * particle_spacing, 0.0) +inlet = RectangularTank(particle_spacing, open_boundary_size, open_boundary_size, + fluid_density; n_layers=boundary_layers, + velocity=prescribed_velocity * flow_direction, + min_coordinates=min_coords_inlet, + faces=(false, false, true, true), + coordinates_eltype=Float64) + +min_coords_outlet = (pipe.fluid_size[1], 0.0) +outlet = RectangularTank(particle_spacing, open_boundary_size, open_boundary_size, + fluid_density; n_layers=boundary_layers, + velocity=prescribed_velocity * flow_direction, + min_coordinates=min_coords_outlet, + faces=(false, false, true, true), + coordinates_eltype=Float64) + +n_buffer_particles = nparticles(inlet.fluid) + +cylinder_center = (5 * cylinder_diameter, domain_size[2] / 2) +cylinder = SphereShape(particle_spacing, cylinder_diameter / 2, + cylinder_center, fluid_density; sphere_type=RoundSphere()) + +fluid = setdiff(pipe.fluid, cylinder) + +# ========================================================================================== +# ==== Fluid +smoothing_length = 1.5 * particle_spacing +smoothing_kernel = WendlandC2Kernel{2}() + +fluid_density_calculator = ContinuityDensity() + +kinematic_viscosity = prescribed_velocity * cylinder_diameter / reynolds_number + +state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=1) +viscosity = ViscosityAdami(nu=kinematic_viscosity) +density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1) + +shifting_technique = ParticleShiftingTechnique(; sound_speed_factor=0.2, v_max_factor=0) + +fluid_system = WeaklyCompressibleSPHSystem(fluid; smoothing_kernel, smoothing_length, + density_calculator=fluid_density_calculator, + state_equation, density_diffusion, viscosity, + pressure_acceleration=tensile_instability_control, + shifting_technique, + buffer_size=n_buffer_particles) + +# ========================================================================================== +# ==== Open Boundary +open_boundary_model = BoundaryModelMirroringTafuni(; mirror_method=ZerothOrderMirroring()) + +face_in = ([0.0, 0.0], [0.0, domain_size[2]]) +inflow = BoundaryZone(; boundary_face=face_in, face_normal=flow_direction, + open_boundary_layers, density=fluid_density, particle_spacing, + reference_velocity=prescribed_velocity * flow_direction, + initial_condition=inlet.fluid) + +face_out = ([min_coords_outlet[1], 0.0], [min_coords_outlet[1], domain_size[2]]) +outflow = BoundaryZone(; boundary_face=face_out, face_normal=(-flow_direction), + open_boundary_layers, density=fluid_density, particle_spacing, + initial_condition=outlet.fluid) + +open_boundary = OpenBoundarySystem(inflow, outflow; fluid_system, + boundary_model=open_boundary_model, + pressure_acceleration=TrixiParticles.inter_particle_averaged_pressure, + buffer_size=n_buffer_particles) + +# ========================================================================================== +# ==== Boundary +wall = union(pipe.boundary, inlet.boundary, outlet.boundary) +# Free-slip boundary condition for the wall. +boundary_model_wall = BoundaryModelDummyParticles(wall.density, wall.mass, + AdamiPressureExtrapolation(), + smoothing_kernel, smoothing_length; + state_equation) + +boundary_system_wall = WallBoundarySystem(wall, boundary_model_wall) + +# No-slip boundary condition for the cylinder. +boundary_model_cylinder = BoundaryModelDummyParticles(cylinder.density, cylinder.mass, + AdamiPressureExtrapolation(), + smoothing_kernel, smoothing_length; + state_equation, viscosity) + +boundary_system_cylinder = WallBoundarySystem(cylinder, boundary_model_cylinder) + +# ========================================================================================== +# ==== Simulation +min_corner = minimum(wall.coordinates .- particle_spacing, dims=2) +max_corner = maximum(wall.coordinates .+ particle_spacing, dims=2) + +nhs = GridNeighborhoodSearch{2}(; cell_list=FullGridCellList(; min_corner, max_corner), + update_strategy=ParallelUpdate()) + +semi = Semidiscretization(fluid_system, open_boundary, boundary_system_wall, + boundary_system_cylinder; neighborhood_search=nhs, + parallelization_backend=PolyesterBackend()) + +ode = semidiscretize(semi, tspan) + +info_callback = InfoCallback(interval=50) +saving_callback = SolutionSavingCallback(dt=0.02, prefix="") + +extra_callback = nothing + +callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback(), extra_callback) + +sol = solve(ode, RDPK3SpFSAL35(), + abstol=1e-6, # May need tuning to prevent boundary penetration + reltol=1e-4, # May need tuning to prevent boundary penetration + save_everystep=false, callback=callbacks); diff --git a/test/examples/examples_fluid.jl b/test/examples/examples_fluid.jl index f393936209..baecabdcb4 100644 --- a/test/examples/examples_fluid.jl +++ b/test/examples/examples_fluid.jl @@ -525,6 +525,15 @@ @test count_rhs_allocations(sol) == 0 end + @trixi_testset "fluid/vortex_street_2d.jl" begin + @trixi_test_nowarn trixi_include(@__MODULE__, + joinpath(examples_dir(), "fluid", + "vortex_street_2d.jl"), + tspan=(0.0, 0.2)) + @test sol.retcode == ReturnCode.Success + @test count_rhs_allocations(sol) == 0 + end + @trixi_testset "fluid/lid_driven_cavity_2d.jl (EDAC)" begin @trixi_test_nowarn trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid",