Poiseuille flow 3d setup

examples/fluid/poiseuille_flow_2d.jl

@@ -15,10 +15,10 @@ using OrdinaryDiffEq
 
 # ==========================================================================================
 # ==== Resolution
-const wall_distance = 0.001 # distance between top and bottom wall
-const flow_length = 0.004 # distance between inflow and outflow
+channel_height = 0.001 # distance between top and bottom walls
+channel_length = 0.004 # distance between inlet and outlet
 
-particle_spacing = wall_distance / 50
+particle_spacing = channel_height / 30
 
 # Make sure that the kernel support of fluid particles at a boundary is always fully sampled
 boundary_layers = 4
@@ -29,51 +29,53 @@ open_boundary_layers = 10
 
 # ==========================================================================================
 # ==== Experiment Setup
-tspan = (0.0, 2.0)
-wcsph = true
+tspan = (0.0, 1.0)
+use_wcsph = true
 
-domain_size = (flow_length, wall_distance)
+domain_size = (channel_length, channel_height)
 
 open_boundary_size = (open_boundary_layers * particle_spacing, domain_size[2])
 
 fluid_density = 1000.0
 reynolds_number = 50
-const pressure_drop = 0.1
-pressure_out = 0.1
-pressure_in = pressure_out + pressure_drop
-const dynamic_viscosity = sqrt(fluid_density * wall_distance^3 * pressure_drop /
-                               (8 * flow_length * reynolds_number))
+imposed_pressure_drop = 0.1
+outlet_pressure = 0.1
+inlet_pressure = outlet_pressure + imposed_pressure_drop
+const dynamic_viscosity = sqrt(fluid_density * channel_height^3 * imposed_pressure_drop /
+                               (8 * channel_length * reynolds_number))
 
-v_max = wall_distance^2 * pressure_drop / (8 * dynamic_viscosity * flow_length)
+v_max = channel_height^2 * imposed_pressure_drop / (8 * dynamic_viscosity * channel_length)
 
 sound_speed_factor = 100
 sound_speed = sound_speed_factor * v_max
 
 flow_direction = (1.0, 0.0)
 
-pipe = RectangularTank(particle_spacing, domain_size, domain_size, fluid_density,
-                       pressure=(pos) -> pressure_out +
-                                         pressure_drop * (1 - (pos[1] / flow_length)),
-                       n_layers=boundary_layers, faces=(false, false, true, true),
-                       coordinates_eltype=Float64)
+channel = RectangularTank(particle_spacing, domain_size, domain_size, fluid_density,
+                          pressure=(pos) -> outlet_pressure +
+                                            imposed_pressure_drop *
+                                            (1 - (pos[1] / channel_length)),
+                          n_layers=boundary_layers, faces=(false, false, true, true),
+                          coordinates_eltype=Float64)
 
 # The analytical solution depends on the length of the fluid domain.
 # Thus, the `BoundaryZone`s extend into the fluid domain because the pressure is not
 # prescribed at the `boundary_face`, but at the free surface of the `BoundaryZone`.
 inlet = RectangularTank(particle_spacing, open_boundary_size, open_boundary_size,
-                        fluid_density, n_layers=boundary_layers, pressure=pressure_in,
+                        fluid_density, n_layers=boundary_layers, pressure=inlet_pressure,
                         min_coordinates=(0.0, 0.0), faces=(false, false, true, true),
                         coordinates_eltype=Float64)
 
 outlet = RectangularTank(particle_spacing, open_boundary_size, open_boundary_size,
                          fluid_density, n_layers=boundary_layers,
-                         min_coordinates=(pipe.fluid_size[1] - open_boundary_size[1], 0.0),
+                         min_coordinates=(channel.fluid_size[1] - open_boundary_size[1],
+                                          0.0),
                          faces=(false, false, true, true),
                          coordinates_eltype=Float64)
 
-fluid = setdiff(pipe.fluid, inlet.fluid, outlet.fluid)
+fluid = setdiff(channel.fluid, inlet.fluid, outlet.fluid)
 
-n_buffer_particles = 10 * pipe.n_particles_per_dimension[2]^2
+n_buffer_particles = 10 * channel.n_particles_per_dimension[2]^2
 
 # ==========================================================================================
 # ==== Fluid
@@ -89,7 +91,7 @@ viscosity = ViscosityAdami(nu=kinematic_viscosity)
 background_pressure = 7 * sound_speed_factor / 10 * fluid_density * v_max^2
 shifting_technique = TransportVelocityAdami(; background_pressure)
 
-if wcsph
+if use_wcsph
     state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density,
                                        exponent=1)
     fluid_system = WeaklyCompressibleSPHSystem(fluid, fluid_density_calculator,
@@ -113,48 +115,54 @@ end
 # ==== Open Boundary
 open_boundary_model = BoundaryModelDynamicalPressureZhang()
 
-boundary_type_in = BidirectionalFlow()
-face_in = ([open_boundary_size[1], 0.0], [open_boundary_size[1], pipe.fluid_size[2]])
-reference_velocity_in = nothing
-reference_pressure_in = 0.2
-inflow = BoundaryZone(; boundary_face=face_in, face_normal=flow_direction,
-                      open_boundary_layers, density=fluid_density, particle_spacing,
-                      reference_velocity=reference_velocity_in,
-                      reference_pressure=reference_pressure_in,
-                      initial_condition=inlet.fluid, boundary_type=boundary_type_in)
-
-boundary_type_out = BidirectionalFlow()
-face_out = ([pipe.fluid_size[1] - open_boundary_size[1], 0.0],
-            [pipe.fluid_size[1] - open_boundary_size[1], pipe.fluid_size[2]])
-reference_velocity_out = nothing
-reference_pressure_out = 0.1
-outflow = BoundaryZone(; boundary_face=face_out, face_normal=(.-(flow_direction)),
-                       open_boundary_layers, density=fluid_density, particle_spacing,
-                       reference_velocity=reference_velocity_out,
-                       reference_pressure=reference_pressure_out,
-                       initial_condition=outlet.fluid, boundary_type=boundary_type_out)
-
-open_boundary = OpenBoundarySystem(inflow, outflow; fluid_system,
+inlet_boundary_type = BidirectionalFlow()
+inlet_face = ([open_boundary_size[1], 0.0], [open_boundary_size[1], channel.fluid_size[2]])
+inlet_reference_velocity = nothing
+inlet_reference_pressure = 0.2
+inlet_boundary_zone = BoundaryZone(; boundary_face=inlet_face, face_normal=flow_direction,
+                                   open_boundary_layers, density=fluid_density,
+                                   particle_spacing,
+                                   reference_velocity=inlet_reference_velocity,
+                                   reference_pressure=inlet_reference_pressure,
+                                   initial_condition=inlet.fluid,
+                                   boundary_type=inlet_boundary_type)
+
+outlet_boundary_type = BidirectionalFlow()
+outlet_face = ([channel.fluid_size[1] - open_boundary_size[1], 0.0],
+               [channel.fluid_size[1] - open_boundary_size[1], channel.fluid_size[2]])
+outlet_reference_velocity = nothing
+outlet_reference_pressure = 0.1
+outlet_flow_direction = (.-(flow_direction))
+outlet_boundary_zone = BoundaryZone(; boundary_face=outlet_face,
+                                    face_normal=outlet_flow_direction,
+                                    open_boundary_layers, density=fluid_density,
+                                    particle_spacing,
+                                    reference_velocity=outlet_reference_velocity,
+                                    reference_pressure=outlet_reference_pressure,
+                                    initial_condition=outlet.fluid,
+                                    boundary_type=outlet_boundary_type)
+
+open_boundary = OpenBoundarySystem(inlet_boundary_zone, outlet_boundary_zone; fluid_system,
                                    boundary_model=open_boundary_model,
                                    calculate_flow_rate=true,
                                    buffer_size=n_buffer_particles)
 
 # ==========================================================================================
 # ==== Boundary
-wall = union(pipe.boundary)
+wall_boundary = union(channel.boundary)
 
-boundary_model = BoundaryModelDummyParticles(wall.density, wall.mass,
+boundary_model = BoundaryModelDummyParticles(wall_boundary.density, wall_boundary.mass,
                                              AdamiPressureExtrapolation(),
                                              state_equation=state_equation,
                                              viscosity=viscosity,
                                              smoothing_kernel, smoothing_length)
 
-boundary_system = WallBoundarySystem(wall, boundary_model)
+boundary_system = WallBoundarySystem(wall_boundary, boundary_model)
 
 # ==========================================================================================
 # ==== Simulation
-min_corner = minimum(wall.coordinates .- 2 * particle_spacing, dims=2)
-max_corner = maximum(wall.coordinates .+ 2 * particle_spacing, dims=2)
+min_corner = minimum(wall_boundary.coordinates .- 2 * particle_spacing, dims=2)
+max_corner = maximum(wall_boundary.coordinates .+ 2 * particle_spacing, dims=2)
 
 nhs = GridNeighborhoodSearch{2}(; cell_list=FullGridCellList(; min_corner, max_corner),
                                 update_strategy=ParallelUpdate())

examples/fluid/poiseuille_flow_3d.jl

@@ -0,0 +1,199 @@
+# ==========================================================================================
+# 3D Hagen-Poiseuille Flow Simulation (Weakly Compressible SPH)
+#
+# Based on:
+#   Zhan, X., et al. "Dynamical pressure boundary condition for weakly compressible smoothed particle hydrodynamics"
+#   Physics of Fluids, Volume 37
+#   https://doi.org/10.1063/5.0254575
+#
+# This example sets up a 3D Hagen-Poiseuille flow simulation in a circular pipe
+# including open boundary conditions.
+# ==========================================================================================
+using TrixiParticles
+using OrdinaryDiffEq
+
+# ==========================================================================================
+# ==== Resolution
+channel_length = 0.004
+channel_radius = 0.0005
+channel_diameter = 2 * channel_radius
+
+particle_spacing_factor = 15
+particle_spacing = channel_diameter / particle_spacing_factor
+
+# 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.
+open_boundary_layers = 10
+
+# ==========================================================================================
+# ==== Experiment Setup
+tspan = (0.0, 1.0)
+
+fluid_density = 1000.0
+reynolds_number = 50
+imposed_pressure_drop = 0.1
+const dynamic_viscosity = sqrt(fluid_density * channel_radius^2 * imposed_pressure_drop /
+                               (4 * reynolds_number))
+
+v_max = channel_radius^2 * imposed_pressure_drop / (4 * dynamic_viscosity * channel_length)
+
+sound_speed_factor = 100
+sound_speed = sound_speed_factor * v_max
+
+flow_direction = (1.0, 0.0, 0.0)
+outlet_flow_direction = (-flow_direction[1], -flow_direction[2], -flow_direction[3])
+
+n_particles_length = ceil(Int, channel_length / particle_spacing)
+
+wall_cross_section = SphereShape(particle_spacing, channel_radius, (0.0, 0.0),
+                                 fluid_density, n_layers=boundary_layers,
+                                 layer_outwards=true, sphere_type=RoundSphere())
+
+# Extend 2D coordinates to 3D by adding x-coordinates
+wall_cross_section_coordinates = hcat(particle_spacing / 2 *
+                                      ones(nparticles(wall_cross_section)),
+                                      wall_cross_section.coordinates')'
+
+wall_boundary = extrude_geometry(wall_cross_section_coordinates; particle_spacing,
+                                 direction=collect(flow_direction),
+                                 density=fluid_density,
+                                 n_extrude=n_particles_length)
+
+fluid_cross_section = SphereShape(particle_spacing, channel_radius, (0.0, 0.0),
+                                  fluid_density, sphere_type=RoundSphere())
+
+# Extend 2D coordinates to 3D by adding x-coordinates
+fluid_cross_section_coordinates = hcat(particle_spacing / 2 *
+                                       ones(nparticles(fluid_cross_section)),
+                                       fluid_cross_section.coordinates')'
+
+inlet_particles = extrude_geometry(fluid_cross_section_coordinates; particle_spacing,
+                                   direction=collect(flow_direction),
+                                   density=fluid_density,
+                                   n_extrude=open_boundary_layers)
+
+fluid_offset = [open_boundary_layers * particle_spacing, 0, 0]
+fluid_particles = extrude_geometry(fluid_cross_section_coordinates .+ fluid_offset;
+                                   particle_spacing,
+                                   direction=collect(flow_direction),
+                                   density=fluid_density,
+                                   n_extrude=(n_particles_length - 2 * open_boundary_layers))
+
+outlet_offset = [(n_particles_length - open_boundary_layers) * particle_spacing, 0, 0]
+outlet_particles = extrude_geometry(fluid_cross_section_coordinates .+ outlet_offset;
+                                    particle_spacing,
+                                    direction=collect(flow_direction),
+                                    n_extrude=open_boundary_layers,
+                                    density=fluid_density)
+
+n_buffer_particles = 10 * nparticles(fluid_cross_section)
+
+# ==========================================================================================
+# ==== Fluid
+smoothing_length = 2 * particle_spacing
+smoothing_kernel = WendlandC2Kernel{3}()
+
+fluid_density_calculator = ContinuityDensity()
+
+kinematic_viscosity = dynamic_viscosity / fluid_density
+viscosity = ViscosityAdami(nu=kinematic_viscosity)
+
+background_pressure = 7 * sound_speed_factor / 10 * fluid_density * v_max^2
+shifting_technique = TransportVelocityAdami(; background_pressure)
+
+state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density,
+                                   exponent=1)
+
+fluid_system = WeaklyCompressibleSPHSystem(fluid_particles, fluid_density_calculator,
+                                           state_equation, smoothing_kernel,
+                                           buffer_size=n_buffer_particles,
+                                           shifting_technique=shifting_technique,
+                                           density_diffusion=DensityDiffusionMolteniColagrossi(delta=0.1),
+                                           smoothing_length, viscosity=viscosity)
+
+# ==========================================================================================
+# ==== Open Boundary
+open_boundary_model = BoundaryModelDynamicalPressureZhang()
+
+inlet_boundary_type = BidirectionalFlow()
+inlet_face = ([
+                  open_boundary_layers * particle_spacing,
+                  -channel_diameter,
+                  -channel_diameter
+              ],
+              [
+                  open_boundary_layers * particle_spacing,
+                  -channel_diameter,
+                  channel_diameter
+              ],
+              [open_boundary_layers * particle_spacing, channel_diameter, channel_diameter])
+inlet_reference_velocity = nothing
+inlet_reference_pressure = 0.2
+
+inlet_zone = BoundaryZone(; boundary_face=inlet_face, face_normal=flow_direction,
+                          open_boundary_layers, density=fluid_density, particle_spacing,
+                          reference_velocity=inlet_reference_velocity,
+                          reference_pressure=inlet_reference_pressure,
+                          initial_condition=inlet_particles,
+                          boundary_type=inlet_boundary_type)
+
+outlet_boundary_type = BidirectionalFlow()
+outlet_face = ([outlet_offset[1], -channel_diameter, -channel_diameter],
+               [outlet_offset[1], -channel_diameter, channel_diameter],
+               [outlet_offset[1], channel_diameter, channel_diameter])
+outlet_reference_velocity = nothing
+outlet_reference_pressure = 0.1
+
+outlet_zone = BoundaryZone(; boundary_face=outlet_face,
+                           face_normal=outlet_flow_direction,
+                           open_boundary_layers, density=fluid_density, particle_spacing,
+                           reference_velocity=outlet_reference_velocity,
+                           reference_pressure=outlet_reference_pressure,
+                           initial_condition=outlet_particles,
+                           boundary_type=outlet_boundary_type)
+
+open_boundary = OpenBoundarySystem(inlet_zone, outlet_zone; fluid_system,
+                                   boundary_model=open_boundary_model,
+                                   buffer_size=n_buffer_particles)
+
+# ==========================================================================================
+# ==== Boundary
+boundary_model = BoundaryModelDummyParticles(wall_boundary.density, wall_boundary.mass,
+                                             AdamiPressureExtrapolation(),
+                                             state_equation=state_equation,
+                                             viscosity=viscosity,
+                                             smoothing_kernel, smoothing_length)
+
+boundary_system = WallBoundarySystem(wall_boundary, boundary_model)
+
+# ==========================================================================================
+# ==== Simulation
+min_corner = minimum(wall_boundary.coordinates .- 2 * particle_spacing, dims=2)
+max_corner = maximum(wall_boundary.coordinates .+ 2 * particle_spacing, dims=2)
+
+neighborhood_search = GridNeighborhoodSearch{3}(;
+                                                cell_list=FullGridCellList(; min_corner,
+                                                                           max_corner),
+                                                update_strategy=ParallelUpdate())
+
+semi_discretization = Semidiscretization(fluid_system, open_boundary,
+                                         boundary_system,
+                                         neighborhood_search=neighborhood_search,
+                                         parallelization_backend=PolyesterBackend())
+
+ode_problem = semidiscretize(semi_discretization, tspan)
+
+info_callback = InfoCallback(interval=20)
+saving_callback = SolutionSavingCallback(dt=0.01, prefix="", output_directory="out")
+extra_callback = nothing
+
+callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback(), extra_callback)
+
+sol = solve(ode_problem, RDPK3SpFSAL35(),
+            abstol=1e-6, # Default abstol is 1e-6 (may need tuning to prevent boundary penetration)
+            reltol=1e-4, # Default reltol is 1e-3 (may need tuning to prevent boundary penetration)
+            dtmax=1e-2,  # Limit stepsize to prevent crashing
+            save_everystep=false, callback=callbacks)