Validate open boundaries

validation/vortex_street_2d/plot_vortex_street_2d.jl

@@ -0,0 +1,2 @@
+# TODO: Plot the comparison with the reference solution
+# C_L and C_D from Tafuni et al. (2018)

validation/vortex_street_2d/validation_vortex_street_2d.jl

@@ -0,0 +1,63 @@
+using TrixiParticles
+using FFTW
+using CSV, DataFrames
+using Test
+
+# Results in [90k particles, 220k particles, 1.2M particles, 5M particles]
+# In the Tafuni et al. (2018), the resolution is `0.01` (5M particles).
+resolution_factor = 0.08 # [0.08, 0.05, 0.02, 0.01]
+
+# ======================================================================================
+# ==== Run the simulation
+trixi_include(joinpath(validation_dir(), "vortex_street_2d", "vortex_street_2d.jl"),
+              parallelization_backend=PolyesterBackend(),
+              factor_d=resolution_factor, saving_callback=nothing, tspan=(0.0, 20.0))
+
+# ======================================================================================
+# ==== Read results and compute the Strouhal number
+data = CSV.read(joinpath(output_directory, "resulting_force.csv"), DataFrame)
+
+time = data[!, "time"]
+t_start = 6.0
+start_index = findfirst(t -> t ≥ t_start, time)
+times_cut = time[start_index:end]
+dt = times_cut[2] - times_cut[1]
+
+f_lift = data[!, "f_l_fluid_1"][start_index:end]
+
+# Compute the frequency for the FFT.
+# For N time samples with uniform time steps dt, the corresponding frequencies are:
+# f_k = k / (N * dt), where k = 0, 1, ..., N-1.
+# This gives the frequency bins in Hz, matching the order of FFT.
+N = length(f_lift)
+frequencies = (0:(N - 1)) / (N * dt)
+
+spectrum = abs.(fft(f_lift))
+
+# For real-valued signals, the FFT output is symmetric.
+# Only the first half (up to the Nyquist frequency) contains unique, physically meaningful frequency components.
+# We therefore analyze only the first N/2 values of the frequency spectrum.
+f_dominant = frequencies[argmax(spectrum[1:div(N, 2)])]
+
+# Verify whether the dominant frequency is indeed unique.
+# In theory, for a purely harmonic oscillation, the spectrum should exhibit only a single dominant frequency component.
+delta = 2 * (frequencies[2] - frequencies[1])
+frequency_band = (abs.(frequencies[1:div(N, 2)] .- f_dominant) .< delta)
+
+# ======================================================================================
+# ==== Save the strouhal numbers
+strouhal_number = f_dominant * cylinder_diameter / prescribed_velocity
+
+df = DataFrame(Resolution=resolution_factor, t_max=last(tspan),
+               frequency=f_dominant, StrouhalNumber=strouhal_number)
+
+CSV.write(joinpath(output_directory, "strouhal_number.csv"), df)
+
+# ======================================================================================
+# ==== Validate the frequency spectrum
+spectrum_half = spectrum[1:div(N, 2)]
+integral_total = sum(spectrum_half)
+integral_peak = sum(spectrum_half[frequency_band])
+
+# TODO: 0.4 is sufficient? Check for higher resolution.
+@test 0.4 < integral_peak / integral_total

validation/vortex_street_2d/vortex_street_2d.jl

@@ -0,0 +1,225 @@
+# Flow past a circular cylinder (vortex street), Tafuni et al. (2018).
+# Other literature using this validation:
+# Vacandio et al. (2013), Marrone et al. (2013), Calhoun (2002), Liu et al. (1998)
+
+using TrixiParticles
+using OrdinaryDiffEq
+
+# ==========================================================================================
+# ==== Resolution
+factor_d = 0.08 # Resolution in the paper is `0.01` (5M particles)
+
+const cylinder_diameter = 0.1
+particle_spacing = factor_d * 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, 20.0)
+
+# Boundary geometry and initial fluid particle positions
+domain_size = (25 * cylinder_diameter, 20 * cylinder_diameter)
+
+flow_direction = [1.0, 0.0]
+reynolds_number = 200
+const prescribed_velocity = 1.0
+const fluid_density = 1000.0
+sound_speed = 10 * prescribed_velocity
+
+boundary_size = (domain_size[1] + 2 * particle_spacing * open_boundary_layers,
+                 domain_size[2])
+
+pipe = RectangularTank(particle_spacing, domain_size, boundary_size, fluid_density,
+                       n_layers=boundary_layers, velocity=[prescribed_velocity, 0.0],
+                       faces=(false, false, true, true))
+
+# Shift pipe walls in negative x-direction for the inflow
+pipe.boundary.coordinates[1, :] .-= particle_spacing * open_boundary_layers
+
+n_buffer_particles = 10 * pipe.n_particles_per_dimension[2]^(ndims(pipe.fluid) - 1)
+
+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
+wcsph = true
+
+h_factor = 2.0
+smoothing_length = h_factor * 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=7)
+viscosity = ViscosityAdami(nu=kinematic_viscosity)
+density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1)
+
+fluid_system = WeaklyCompressibleSPHSystem(fluid, fluid_density_calculator,
+                                           state_equation, smoothing_kernel,
+                                           density_diffusion=density_diffusion,
+                                           smoothing_length, viscosity=viscosity,
+                                           pressure_acceleration=tensile_instability_control,
+                                           buffer_size=n_buffer_particles)
+
+# ==========================================================================================
+# ==== Open Boundary
+function velocity_function2d(pos, t)
+    return SVector(prescribed_velocity, 0.0)
+end
+
+open_boundary_model = BoundaryModelTafuni()
+
+boundary_type_in = InFlow()
+plane_in = ([0.0, 0.0], [0.0, domain_size[2]])
+inflow = BoundaryZone(; plane=plane_in, plane_normal=flow_direction, open_boundary_layers,
+                      density=fluid_density, particle_spacing,
+                      boundary_type=boundary_type_in)
+
+reference_velocity_in = velocity_function2d
+# At the inlet, neither pressure nor density are prescribed; instead,
+# these values are extrapolated from the fluid domain
+reference_pressure_in = nothing
+reference_density_in = nothing
+open_boundary_in = OpenBoundarySPHSystem(inflow; fluid_system,
+                                         boundary_model=open_boundary_model,
+                                         buffer_size=n_buffer_particles,
+                                         reference_density=reference_density_in,
+                                         reference_pressure=reference_pressure_in,
+                                         reference_velocity=reference_velocity_in)
+
+boundary_type_out = OutFlow()
+plane_out = ([domain_size[1], 0.0], [domain_size[1], domain_size[2]])
+outflow = BoundaryZone(; plane=plane_out, plane_normal=(-flow_direction),
+                       open_boundary_layers, density=fluid_density, particle_spacing,
+                       boundary_type=boundary_type_out)
+
+# At the outlet, we allow the flow to exit freely without imposing any boundary conditions
+reference_velocity_out = nothing
+reference_pressure_out = nothing
+reference_density_out = nothing
+open_boundary_out = OpenBoundarySPHSystem(outflow; fluid_system,
+                                          boundary_model=open_boundary_model,
+                                          buffer_size=n_buffer_particles,
+                                          reference_density=reference_density_out,
+                                          reference_pressure=reference_pressure_out,
+                                          reference_velocity=reference_velocity_out)
+# ==========================================================================================
+# ==== Boundary
+boundary_model = BoundaryModelDummyParticles(pipe.boundary.density, pipe.boundary.mass,
+                                             AdamiPressureExtrapolation(),
+                                             state_equation=state_equation,
+                                             smoothing_kernel, smoothing_length)
+
+boundary_system_wall = BoundarySPHSystem(pipe.boundary, boundary_model)
+
+boundary_model_cylinder = BoundaryModelDummyParticles(cylinder.density, cylinder.mass,
+                                                      AdamiPressureExtrapolation(),
+                                                      state_equation=state_equation,
+                                                      viscosity=viscosity,
+                                                      smoothing_kernel, smoothing_length)
+
+boundary_system_cylinder = BoundarySPHSystem(cylinder, boundary_model_cylinder)
+
+# ==========================================================================================
+# ==== Postprocessing
+circle = SphereShape(particle_spacing, (cylinder_diameter + particle_spacing) / 2,
+                     cylinder_center, fluid_density, n_layers=1,
+                     sphere_type=RoundSphere())
+
+# Points for pressure interpolation, located at the wall interface
+const data_points =  copy(circle.coordinates)
+const center = SVector(cylinder_center)
+
+calculate_lift_force(system, v_ode, u_ode, semi, t) = nothing
+function calculate_lift_force(system::TrixiParticles.FluidSystem, v_ode, u_ode, semi, t)
+    force = zero(SVector{ndims(system), eltype(system)})
+
+    values = interpolate_points(data_points, semi, system, v_ode, u_ode; cut_off_bnd=false,
+                                clip_negative_pressure=false)
+    pressure = Array(values.pressure)
+
+    for i in axes(data_points, 2)
+        point = TrixiParticles.current_coords(data_points, system, i)
+
+        # F = ∑ -p_i * A_i * n_i
+        force -= pressure[i] * particle_spacing .*
+                 TrixiParticles.normalize(point - center)
+    end
+
+    return 2 * force[2] / (fluid_density * prescribed_velocity^2 * cylinder_diameter)
+end
+
+calculate_drag_force(system, v_ode, u_ode, semi, t) = nothing
+function calculate_drag_force(system::TrixiParticles.FluidSystem, v_ode, u_ode, semi, t)
+    force = zero(SVector{ndims(system), eltype(system)})
+
+    values = interpolate_points(data_points, semi, system, v_ode, u_ode; cut_off_bnd=false,
+                                clip_negative_pressure=false)
+    pressure = Array(values.pressure)
+
+    for i in axes(data_points, 2)
+        point = TrixiParticles.current_coords(data_points, system, i)
+
+        # F = ∑ -p_i * A_i * n_i
+        force -= pressure[i] * particle_spacing .*
+                 TrixiParticles.normalize(point - center)
+    end
+
+    return 2 * force[1] / (fluid_density * prescribed_velocity^2 * cylinder_diameter)
+end
+
+# ==========================================================================================
+# ==== Simulation
+min_corner = minimum(pipe.boundary.coordinates .- particle_spacing, dims=2)
+max_corner = maximum(pipe.boundary.coordinates .+ particle_spacing, dims=2)
+cell_list = FullGridCellList(; min_corner, max_corner)
+
+neighborhood_search = GridNeighborhoodSearch{2}(; cell_list,
+                                                update_strategy=ParallelUpdate())
+
+semi = Semidiscretization(fluid_system, open_boundary_in, open_boundary_out,
+                          boundary_system_wall, boundary_system_cylinder;
+                          neighborhood_search=neighborhood_search,
+                          parallelization_backend=PolyesterBackend())
+
+ode = semidiscretize(semi, tspan)
+
+info_callback = InfoCallback(interval=100)
+
+output_directory = joinpath(validation_dir(), "vortex_street_2d",
+                            "out_vortex_street_dp_$(factor_d)D_c_$(sound_speed)_h_factor_$(h_factor)_" *
+                            TrixiParticles.type2string(smoothing_kernel))
+
+saving_callback = SolutionSavingCallback(; dt=0.02, prefix="", output_directory)
+
+pp_callback = PostprocessCallback(; dt=0.02,
+                                  f_l=calculate_lift_force, f_d=calculate_drag_force,
+                                  output_directory, filename="resulting_force",
+                                  write_csv=true, write_file_interval=10)
+
+extra_callback = nothing
+
+callbacks = CallbackSet(info_callback, UpdateCallback(), saving_callback,
+                        ParticleShiftingCallback(), # To obtain a near-uniform particle distribution in the wake
+                        pp_callback, extra_callback)
+
+sol = solve(ode, RDPK3SpFSAL35(),
+            abstol=1e-6, # 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);