Skip to content

Commit

Permalink
hopefully final fix
Browse files Browse the repository at this point in the history 
The `reference_velocity` fixes the issue with the "wrong velocity" when
domain velocity is not equal the inflow velocity
  • Loading branch information
@LasNikas
LasNikas committed Dec 6, 2023
1 parent 71e2ee5 commit feaba99
Show file tree
Hide file tree
Showing 8 changed files with 110 additions and 234 deletions.
47 changes: 30 additions & 17 deletions examples/fluid/bent_pipe_2d.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,8 @@
using TrixiParticles
using OrdinaryDiffEq

wcsph = false

# ==========================================================================================
# ==== Resolution
particle_spacing = 0.05
Expand All @@ -23,10 +25,10 @@ prescribed_velocity = (0.0, 2.0)

reynolds_number = 1000.0
fluid_density = 1000.0
pressure = 0.0
pressure = wcsph ? 10_000.0 : 100_000.0
sound_speed = 10 * maximum(prescribed_velocity)

pipe_size = (pipe_radius - pipe_radius_inner, 2pipe_radius)
pipe_size = (pipe_radius - pipe_radius_inner, pipe_radius)
pipe_coords_min = (-(boundary_layers * particle_spacing + pipe_radius), 0.0)
pipe_coords_max = (boundary_layers * particle_spacing + 2pipe_radius, pipe_size[2])

Expand Down Expand Up @@ -58,16 +60,24 @@ n_buffer_particles = 20 * pipe_in.n_particles_per_dimension[1]
smoothing_length = 1.2 * particle_spacing
smoothing_kernel = SchoenbergQuinticSplineKernel{2}()

fluid_density_calculator = ContinuityDensity()
viscosity = ArtificialViscosityMonaghan(alpha=0.02, beta=0.0)
state_equation = StateEquationCole(sound_speed, 7, fluid_density, pressure,
background_pressure=pressure)

fluid_system = WeaklyCompressibleSPHSystem(fluid, fluid_density_calculator,
state_equation, smoothing_kernel,
smoothing_length, viscosity=viscosity,
buffer=n_buffer_particles)

if wcsph
fluid_density_calculator = ContinuityDensity()
viscosity = ArtificialViscosityMonaghan(alpha=0.02, beta=0.0)
state_equation = StateEquationCole(sound_speed, 7, fluid_density, pressure)

fluid_system = WeaklyCompressibleSPHSystem(fluid, fluid_density_calculator,
state_equation, smoothing_kernel,
smoothing_length, viscosity=viscosity,
buffer=n_buffer_particles)
else
nu = maximum(prescribed_velocity) * (pipe_radius - pipe_radius_inner) / reynolds_number
viscosity = ViscosityAdami(; nu) #alpha * smoothing_length * sound_speed / 8)

fluid_system = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, smoothing_length,
sound_speed, viscosity=viscosity,
transport_velocity=TransportVelocityAdami(pressure),
buffer=n_buffer_particles)
end
# ==========================================================================================
# ==== Open Boundary
open_boundary_length = open_boundary_cols * particle_spacing
Expand All @@ -91,21 +101,21 @@ zone_plane_in = ([0.0; 0.0], [pipe_size[1]; 0.0])
zone_plane_out = ([pipe_radius + pipe_radius_inner; -open_boundary_length],
[pipe_radius + pipe_radius_inner + pipe_size[1]; -open_boundary_length])

v_x(position) = prescribed_velocity[1]
v_y(position) = prescribed_velocity[2]

open_boundary_in = OpenBoundarySPHSystem(inflow.fluid, InFlow(), fluid_system,
flow_direction=(0.0, 1.0),
zone_width=open_boundary_length,
zone_plane_min_corner=[0.0, 0.0],
zone_plane_max_corner=[pipe_size[1], 0.0],
buffer=n_buffer_particles,
initial_velocity_function=(v_x, v_y))
buffer=n_buffer_particles)

zone_plane_min_corner = [pipe_radius + pipe_radius_inner, 0.0]
zone_plane_max_corner = [pipe_radius + pipe_radius_inner + pipe_size[1], 0.0]
v_x(position, t) = prescribed_velocity[1]
v_y_out(position, t) = -prescribed_velocity[2]

open_boundary_out = OpenBoundarySPHSystem(outflow.fluid, OutFlow(), fluid_system,
flow_direction=(0.0, -1.0),
velocity_function=(v_x, v_y_out),
zone_width=open_boundary_length,
zone_plane_min_corner=zone_plane_min_corner,
zone_plane_max_corner=zone_plane_max_corner,
Expand All @@ -115,9 +125,12 @@ open_boundary_out = OpenBoundarySPHSystem(outflow.fluid, OutFlow(), fluid_system
# ==== Boundary
boundary = union(pipe, inflow.boundary, outflow.boundary)
boundary_density_calculator = AdamiPressureExtrapolation()
state_equation = wcsph ? state_equation : nothing

boundary_model = BoundaryModelDummyParticles(boundary.density, boundary.mass,
boundary_density_calculator,
state_equation=state_equation,
#viscosity=viscosity,
smoothing_kernel, smoothing_length)

boundary_system = BoundarySPHSystem(boundary, boundary_model)
Expand Down
150 changes: 0 additions & 150 deletions examples/fluid/bent_pipe_edac_2d.jl
View file

This file was deleted.

48 changes: 34 additions & 14 deletions examples/fluid/pipe_flow_2d.jl
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
# TODO: Description
using TrixiParticles
using OrdinaryDiffEq
wcsph = false

# ==========================================================================================
# ==== Resolution
Expand All @@ -24,17 +25,18 @@ reynolds_number = 100
particle_spacing = domain_length_factor * domain_length

fluid_density = 1000.0
pressure = 1000.0
pressure = wcsph ? 10_000.0 : 100_000.0

prescribed_velocity = (2.0, 0.0)
prescribed_velocity = (4.0, 0.0)

sound_speed = 10 * maximum(prescribed_velocity)

pipe = RectangularTank(particle_spacing, (domain_length, domain_width),
(domain_length + particle_spacing * open_boundary_cols,
domain_width), fluid_density, pressure=pressure,
n_layers=3, spacing_ratio=1, faces=(false, false, true, true),
init_velocity=prescribed_velocity)
loop_order_fluid=:x_first,
#init_velocity=prescribed_velocity,
n_layers=3, spacing_ratio=1, faces=(false, false, true, true))

n_particles_y = Int(floor(domain_width / particle_spacing))
n_buffer_particles = 4 * n_particles_y
Expand All @@ -44,14 +46,25 @@ n_buffer_particles = 4 * n_particles_y
smoothing_length = 1.2 * particle_spacing
smoothing_kernel = SchoenbergCubicSplineKernel{2}()

nu = maximum(prescribed_velocity) * domain_length / reynolds_number
viscosity = ViscosityAdami(; nu) #alpha * smoothing_length * sound_speed / 8)

fluid_system = EntropicallyDampedSPHSystem(pipe.fluid, smoothing_kernel, smoothing_length,
sound_speed, viscosity=viscosity,
transport_velocity=TransportVelocityAdami(pressure),
buffer=n_buffer_particles)

if wcsph
fluid_density_calculator = ContinuityDensity()
viscosity = ArtificialViscosityMonaghan(alpha=0.02, beta=0.0)
state_equation = StateEquationCole(sound_speed, 7, fluid_density, pressure)

fluid_system = WeaklyCompressibleSPHSystem(pipe.fluid, fluid_density_calculator,
state_equation, smoothing_kernel,
smoothing_length, viscosity=viscosity,
buffer=n_buffer_particles)
else
nu = maximum(prescribed_velocity) * domain_length / reynolds_number
viscosity = ViscosityAdami(; nu) #alpha * smoothing_length * sound_speed / 8)

fluid_system = EntropicallyDampedSPHSystem(pipe.fluid, smoothing_kernel,
smoothing_length,
sound_speed, viscosity=viscosity,
transport_velocity=TransportVelocityAdami(pressure),
buffer=n_buffer_particles)
end
# ==========================================================================================
# ==== Open Boundary
open_boundary_length = particle_spacing * open_boundary_cols
Expand All @@ -62,10 +75,11 @@ inflow = RectangularTank(particle_spacing, open_boundary_size, open_boundary_siz
init_velocity=prescribed_velocity, pressure=pressure,
min_coordinates=(-open_boundary_length, 0.0),
spacing_ratio=spacing_ratio,
loop_order_fluid=:x_first,
faces=(false, false, true, true))
outflow = RectangularTank(particle_spacing, open_boundary_size, open_boundary_size,
fluid_density; n_layers=boundary_layers,
init_velocity=prescribed_velocity, pressure=pressure,
loop_order_fluid=:x_first, pressure=pressure,
min_coordinates=(domain_length, 0.0), spacing_ratio=spacing_ratio,
faces=(false, false, true, true))

Expand All @@ -76,8 +90,12 @@ open_boundary_in = OpenBoundarySPHSystem(inflow.fluid, InFlow(), fluid_system,
zone_plane_max_corner=[0.0, domain_width],
buffer=n_buffer_particles)

v_x(p, t) = prescribed_velocity[1]
v_y(p, t) = prescribed_velocity[2]

open_boundary_out = OpenBoundarySPHSystem(outflow.fluid, OutFlow(), fluid_system,
flow_direction=(1.0, 0.0),
velocity_function=(v_x, v_y),
zone_width=open_boundary_length,
zone_plane_min_corner=[domain_length, 0.0],
zone_plane_max_corner=[domain_length,
Expand All @@ -88,9 +106,11 @@ open_boundary_out = OpenBoundarySPHSystem(outflow.fluid, OutFlow(), fluid_system
# ==== Boundary
boundary = union(pipe.boundary, inflow.boundary, outflow.boundary)

state_equation = wcsph ? state_equation : nothing
boundary_model = BoundaryModelDummyParticles(boundary.density, boundary.mass,
AdamiPressureExtrapolation(),
viscosity=ViscosityAdami(nu=1e-4),
state_equation=state_equation,
#viscosity=ViscosityAdami(nu=1e-4),
smoothing_kernel, smoothing_length)

boundary_system = BoundarySPHSystem(boundary, boundary_model)
Expand Down
5 changes: 2 additions & 3 deletions examples/fluid/taylor_green_vortex_2d.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -133,7 +133,7 @@ saving_callback = SolutionSavingCallback(dt=0.02,
p_avg=average_pressure,
t=time_vector)

callbacks = CallbackSet(info_callback, saving_callback)
callbacks = CallbackSet(info_callback, saving_callback, UpdateEachTimeStep())

# Use a Runge-Kutta method with automatic (error based) time step size control.
# Enable threading of the RK method for better performance on multiple threads.
Expand All @@ -143,8 +143,7 @@ callbacks = CallbackSet(info_callback, saving_callback)
# Sometimes, the method fails to do so with Monaghan-Kajtar BC because forces
# become extremely large when fluid particles are very close to boundary particles,
# and the time integration method interprets this as an instability.
sol = solve(ode,
RDPK3SpFSAL49((step_limiter!)=TrixiParticles.update_transport_velocity!),
sol = solve(ode, RDPK3SpFSAL49(),
abstol=1e-8, # 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=dt_max,#1e-2, # Limit stepsize to prevent crashing
Expand Down
Loading

0 comments on commit feaba99

Please sign in to comment.