diff --git a/validation/vortex_street_2d/plot_vortex_street_2d.jl b/validation/vortex_street_2d/plot_vortex_street_2d.jl new file mode 100644 index 000000000..761b78414 --- /dev/null +++ b/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) diff --git a/validation/vortex_street_2d/validation_vortex_street_2d.jl b/validation/vortex_street_2d/validation_vortex_street_2d.jl new file mode 100644 index 000000000..12585ba0d --- /dev/null +++ b/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 diff --git a/validation/vortex_street_2d/vortex_street_2d.jl b/validation/vortex_street_2d/vortex_street_2d.jl new file mode 100644 index 000000000..b15e8f031 --- /dev/null +++ b/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);