|
| 1 | +# ========================================================================================== |
| 2 | +# 2D Poiseuille Flow Simulation (Weakly Compressible SPH) |
| 3 | +# |
| 4 | +# Based on: |
| 5 | +# Zhan, X., et al. "Dynamical pressure boundary condition for weakly compressible smoothed particle hydrodynamics" |
| 6 | +# Physics of Fluids, Volume 37 |
| 7 | +# https://doi.org/10.1063/5.0254575 |
| 8 | +# |
| 9 | +# This example sets up a 2D Poiseuille flow simulation in a rectangular channel |
| 10 | +# including open boundary conditions. |
| 11 | +# ========================================================================================== |
| 12 | + |
| 13 | +using TrixiParticles |
| 14 | +using OrdinaryDiffEq |
| 15 | + |
| 16 | +# ========================================================================================== |
| 17 | +# ==== Resolution |
| 18 | +const wall_distance = 0.001 # distance between top and bottom wall |
| 19 | +const flow_length = 0.004 # distance between inflow and outflow |
| 20 | + |
| 21 | +particle_spacing = wall_distance / 50 |
| 22 | + |
| 23 | +# Make sure that the kernel support of fluid particles at a boundary is always fully sampled |
| 24 | +boundary_layers = 4 |
| 25 | + |
| 26 | +# Make sure that the kernel support of fluid particles at an open boundary is always |
| 27 | +# fully sampled. |
| 28 | +open_boundary_layers = 10 |
| 29 | + |
| 30 | +# ========================================================================================== |
| 31 | +# ==== Experiment Setup |
| 32 | +tspan = (0.0, 2.0) |
| 33 | +wcsph = true |
| 34 | + |
| 35 | +domain_size = (flow_length, wall_distance) |
| 36 | + |
| 37 | +open_boundary_size = (open_boundary_layers * particle_spacing, domain_size[2]) |
| 38 | + |
| 39 | +fluid_density = 1000.0 |
| 40 | +reynolds_number = 50 |
| 41 | +const pressure_drop = 0.1 |
| 42 | +pressure_out = 0.1 |
| 43 | +pressure_in = pressure_out + pressure_drop |
| 44 | +const dynamic_viscosity = sqrt(fluid_density * wall_distance^3 * pressure_drop / |
| 45 | + (8 * flow_length * reynolds_number)) |
| 46 | + |
| 47 | +v_max = wall_distance^2 * pressure_drop / (8 * dynamic_viscosity * flow_length) |
| 48 | + |
| 49 | +sound_speed_factor = 100 |
| 50 | +sound_speed = sound_speed_factor * v_max |
| 51 | + |
| 52 | +flow_direction = (1.0, 0.0) |
| 53 | + |
| 54 | +pipe = RectangularTank(particle_spacing, domain_size, domain_size, fluid_density, |
| 55 | + pressure=(pos) -> pressure_out + |
| 56 | + pressure_drop * (1 - (pos[1] / flow_length)), |
| 57 | + n_layers=boundary_layers, faces=(false, false, true, true)) |
| 58 | + |
| 59 | +# The analytical solution depends on the length of the fluid domain. |
| 60 | +# Thus, the `BoundaryZone`s extend into the fluid domain because the pressure is not |
| 61 | +# prescribed at the `boundary_face`, but at the free surface of the `BoundaryZone`. |
| 62 | +inlet = RectangularTank(particle_spacing, open_boundary_size, open_boundary_size, |
| 63 | + fluid_density, n_layers=boundary_layers, pressure=pressure_in, |
| 64 | + min_coordinates=(0.0, 0.0), faces=(false, false, true, true)) |
| 65 | + |
| 66 | +outlet = RectangularTank(particle_spacing, open_boundary_size, open_boundary_size, |
| 67 | + fluid_density, n_layers=boundary_layers, |
| 68 | + min_coordinates=(pipe.fluid_size[1] - open_boundary_size[1], 0.0), |
| 69 | + faces=(false, false, true, true)) |
| 70 | + |
| 71 | +fluid = setdiff(pipe.fluid, inlet.fluid, outlet.fluid) |
| 72 | + |
| 73 | +n_buffer_particles = 10 * pipe.n_particles_per_dimension[2]^2 |
| 74 | + |
| 75 | +# ========================================================================================== |
| 76 | +# ==== Fluid |
| 77 | +smoothing_length = 2 * particle_spacing |
| 78 | +smoothing_kernel = WendlandC2Kernel{2}() |
| 79 | + |
| 80 | +fluid_density_calculator = ContinuityDensity() |
| 81 | + |
| 82 | +kinematic_viscosity = dynamic_viscosity / fluid_density |
| 83 | + |
| 84 | +viscosity = ViscosityAdami(nu=kinematic_viscosity) |
| 85 | + |
| 86 | +background_pressure = 7 * sound_speed_factor / 10 * fluid_density * v_max^2 |
| 87 | +shifting_technique = TransportVelocityAdami(; background_pressure) |
| 88 | + |
| 89 | +if wcsph |
| 90 | + state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, |
| 91 | + exponent=1) |
| 92 | + fluid_system = WeaklyCompressibleSPHSystem(fluid, fluid_density_calculator, |
| 93 | + state_equation, smoothing_kernel, |
| 94 | + buffer_size=n_buffer_particles, |
| 95 | + shifting_technique=shifting_technique, |
| 96 | + density_diffusion=DensityDiffusionMolteniColagrossi(delta=0.1), |
| 97 | + smoothing_length, viscosity=viscosity) |
| 98 | +else |
| 99 | + state_equation = nothing |
| 100 | + |
| 101 | + fluid_system = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, |
| 102 | + smoothing_length, |
| 103 | + sound_speed, viscosity=viscosity, |
| 104 | + density_calculator=fluid_density_calculator, |
| 105 | + shifting_technique=shifting_technique, |
| 106 | + buffer_size=n_buffer_particles) |
| 107 | +end |
| 108 | + |
| 109 | +# ========================================================================================== |
| 110 | +# ==== Open Boundary |
| 111 | +open_boundary_model = BoundaryModelDynamicalPressureZhang() |
| 112 | + |
| 113 | +boundary_type_in = BidirectionalFlow() |
| 114 | +face_in = ([open_boundary_size[1], 0.0], [open_boundary_size[1], pipe.fluid_size[2]]) |
| 115 | +reference_velocity_in = nothing |
| 116 | +reference_pressure_in = 0.2 |
| 117 | +inflow = BoundaryZone(; boundary_face=face_in, face_normal=flow_direction, |
| 118 | + open_boundary_layers, density=fluid_density, particle_spacing, |
| 119 | + reference_velocity=reference_velocity_in, |
| 120 | + reference_pressure=reference_pressure_in, |
| 121 | + initial_condition=inlet.fluid, boundary_type=boundary_type_in) |
| 122 | + |
| 123 | +boundary_type_out = BidirectionalFlow() |
| 124 | +face_out = ([pipe.fluid_size[1] - open_boundary_size[1], 0.0], |
| 125 | + [pipe.fluid_size[1] - open_boundary_size[1], pipe.fluid_size[2]]) |
| 126 | +reference_velocity_out = nothing |
| 127 | +reference_pressure_out = 0.1 |
| 128 | +outflow = BoundaryZone(; boundary_face=face_out, face_normal=(.-(flow_direction)), |
| 129 | + open_boundary_layers, density=fluid_density, particle_spacing, |
| 130 | + reference_velocity=reference_velocity_out, |
| 131 | + reference_pressure=reference_pressure_out, |
| 132 | + initial_condition=outlet.fluid, boundary_type=boundary_type_out) |
| 133 | + |
| 134 | +open_boundary = OpenBoundarySystem(inflow, outflow; fluid_system, |
| 135 | + boundary_model=open_boundary_model, |
| 136 | + buffer_size=n_buffer_particles) |
| 137 | + |
| 138 | +# ========================================================================================== |
| 139 | +# ==== Boundary |
| 140 | +wall = union(pipe.boundary) |
| 141 | + |
| 142 | +boundary_model = BoundaryModelDummyParticles(wall.density, wall.mass, |
| 143 | + AdamiPressureExtrapolation(), |
| 144 | + state_equation=state_equation, |
| 145 | + viscosity=viscosity, |
| 146 | + smoothing_kernel, smoothing_length) |
| 147 | + |
| 148 | +boundary_system = WallBoundarySystem(wall, boundary_model) |
| 149 | + |
| 150 | +# ========================================================================================== |
| 151 | +# ==== Simulation |
| 152 | +min_corner = minimum(wall.coordinates .- 2 * particle_spacing, dims=2) |
| 153 | +max_corner = maximum(wall.coordinates .+ 2 * particle_spacing, dims=2) |
| 154 | + |
| 155 | +nhs = GridNeighborhoodSearch{2}(; cell_list=FullGridCellList(; min_corner, max_corner), |
| 156 | + update_strategy=ParallelUpdate()) |
| 157 | + |
| 158 | +semi = Semidiscretization(fluid_system, open_boundary, |
| 159 | + boundary_system, neighborhood_search=nhs, |
| 160 | + parallelization_backend=PolyesterBackend()) |
| 161 | + |
| 162 | +ode = semidiscretize(semi, tspan) |
| 163 | + |
| 164 | +info_callback = InfoCallback(interval=100) |
| 165 | +saving_callback = SolutionSavingCallback(dt=0.02, prefix="", output_directory="out") |
| 166 | + |
| 167 | +extra_callback = nothing |
| 168 | + |
| 169 | +callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback(), extra_callback) |
| 170 | + |
| 171 | +sol = solve(ode, RDPK3SpFSAL35(), |
| 172 | + abstol=1e-6, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration) |
| 173 | + reltol=1e-4, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration) |
| 174 | + dtmax=1e-2, # Limit stepsize to prevent crashing |
| 175 | + save_everystep=false, callback=callbacks); |
0 commit comments