Skip to content

Commit 70fa235

Browse files
svchbsvchb
authored
Poiseuille flow 3d setup (#1075)
* add aqua * Add 3D Hagen-Poiseuille flow simulation example * Add 3D pulsative channel flow simulation example * revert * Refactor Poiseuille flow examples to use consistent variable naming and improve simulation parameters * format * add examples to test * format * shorten sim time * format * replace sol with solution * suppress warnings * fix * format * review * format * fix names * format * use let * add more comments * fix
1 parent 425a7b4 commit 70fa235

File tree

5 files changed

+348
-57
lines changed

5 files changed

+348
-57
lines changed

‎examples/fluid/poiseuille_flow_2d.jl‎

Lines changed: 63 additions & 53 deletions
Original file line numberDiff line numberDiff line change
@@ -15,10 +15,10 @@ using OrdinaryDiffEq
1515

1616
# ==========================================================================================
1717
# ==== Resolution
18-
const wall_distance = 0.001 # distance between top and bottom wall
19-
const flow_length = 0.004 # distance between inflow and outflow
18+
channel_height = 0.001 # distance between top and bottom walls
19+
channel_length = 0.004 # distance between inlet and outlet
2020

21-
particle_spacing = wall_distance / 50
21+
particle_spacing = channel_height / 30
2222

2323
# Make sure that the kernel support of fluid particles at a boundary is always fully sampled
2424
boundary_layers = 4
@@ -29,51 +29,55 @@ open_boundary_layers = 10
2929

3030
# ==========================================================================================
3131
# ==== Experiment Setup
32-
tspan = (0.0, 2.0)
33-
wcsph = true
32+
tspan = (0.0, 1.0)
33+
use_wcsph = true
3434

35-
domain_size = (flow_length, wall_distance)
35+
domain_size = (channel_length, channel_height)
3636

3737
open_boundary_size = (open_boundary_layers * particle_spacing, domain_size[2])
3838

3939
fluid_density = 1000.0
4040
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))
41+
imposed_pressure_drop = 0.1
42+
outlet_pressure = 0.1
43+
inlet_pressure = outlet_pressure + imposed_pressure_drop
44+
const dynamic_viscosity = sqrt(fluid_density * channel_height^3 * imposed_pressure_drop /
45+
(8 * channel_length * reynolds_number))
4646

47-
v_max = wall_distance^2 * pressure_drop / (8 * dynamic_viscosity * flow_length)
47+
v_max = channel_height^2 * imposed_pressure_drop / (8 * dynamic_viscosity * channel_length)
4848

4949
sound_speed_factor = 100
5050
sound_speed = sound_speed_factor * v_max
5151

5252
flow_direction = (1.0, 0.0)
5353

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-
coordinates_eltype=Float64)
59-
60-
# The analytical solution depends on the length of the fluid domain.
61-
# Thus, the `BoundaryZone`s extend into the fluid domain because the pressure is not
62-
# prescribed at the `boundary_face`, but at the free surface of the `BoundaryZone`.
54+
channel = RectangularTank(particle_spacing, domain_size, domain_size, fluid_density,
55+
pressure=(pos) -> outlet_pressure +
56+
imposed_pressure_drop *
57+
(1 - (pos[1] / channel_length)),
58+
n_layers=boundary_layers, faces=(false, false, true, true),
59+
coordinates_eltype=Float64)
60+
61+
# The analytical solution uses the full inflow-to-outflow extent x in [0, channel_length].
62+
# Since the pressure is prescribed at the free surface of each `BoundaryZone` instead of at
63+
# its `boundary_face`, the inlet and outlet free surfaces lie at x = 0 and
64+
# x = channel_length (= L), respectively. Thus, the boundary zones are effectively part of
65+
# the computational domain in this setup.
6366
inlet = RectangularTank(particle_spacing, open_boundary_size, open_boundary_size,
64-
fluid_density, n_layers=boundary_layers, pressure=pressure_in,
67+
fluid_density, n_layers=boundary_layers, pressure=inlet_pressure,
6568
min_coordinates=(0.0, 0.0), faces=(false, false, true, true),
6669
coordinates_eltype=Float64)
6770

6871
outlet = RectangularTank(particle_spacing, open_boundary_size, open_boundary_size,
6972
fluid_density, n_layers=boundary_layers,
70-
min_coordinates=(pipe.fluid_size[1] - open_boundary_size[1], 0.0),
73+
min_coordinates=(channel.fluid_size[1] - open_boundary_size[1],
74+
0.0),
7175
faces=(false, false, true, true),
7276
coordinates_eltype=Float64)
7377

74-
fluid = setdiff(pipe.fluid, inlet.fluid, outlet.fluid)
78+
fluid = setdiff(channel.fluid, inlet.fluid, outlet.fluid)
7579

76-
n_buffer_particles = 10 * pipe.n_particles_per_dimension[2]^2
80+
n_buffer_particles = 10 * channel.n_particles_per_dimension[2]^2
7781

7882
# ==========================================================================================
7983
# ==== Fluid
@@ -89,7 +93,7 @@ viscosity = ViscosityAdami(nu=kinematic_viscosity)
8993
background_pressure = 7 * sound_speed_factor / 10 * fluid_density * v_max^2
9094
shifting_technique = TransportVelocityAdami(; background_pressure)
9195

92-
if wcsph
96+
if use_wcsph
9397
state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density,
9498
exponent=1)
9599
fluid_system = WeaklyCompressibleSPHSystem(fluid, fluid_density_calculator,
@@ -113,48 +117,54 @@ end
113117
# ==== Open Boundary
114118
open_boundary_model = BoundaryModelDynamicalPressureZhang()
115119

116-
boundary_type_in = BidirectionalFlow()
117-
face_in = ([open_boundary_size[1], 0.0], [open_boundary_size[1], pipe.fluid_size[2]])
118-
reference_velocity_in = nothing
119-
reference_pressure_in = 0.2
120-
inflow = BoundaryZone(; boundary_face=face_in, face_normal=flow_direction,
121-
open_boundary_layers, density=fluid_density, particle_spacing,
122-
reference_velocity=reference_velocity_in,
123-
reference_pressure=reference_pressure_in,
124-
initial_condition=inlet.fluid, boundary_type=boundary_type_in)
125-
126-
boundary_type_out = BidirectionalFlow()
127-
face_out = ([pipe.fluid_size[1] - open_boundary_size[1], 0.0],
128-
[pipe.fluid_size[1] - open_boundary_size[1], pipe.fluid_size[2]])
129-
reference_velocity_out = nothing
130-
reference_pressure_out = 0.1
131-
outflow = BoundaryZone(; boundary_face=face_out, face_normal=(.-(flow_direction)),
132-
open_boundary_layers, density=fluid_density, particle_spacing,
133-
reference_velocity=reference_velocity_out,
134-
reference_pressure=reference_pressure_out,
135-
initial_condition=outlet.fluid, boundary_type=boundary_type_out)
136-
137-
open_boundary = OpenBoundarySystem(inflow, outflow; fluid_system,
120+
inlet_boundary_type = BidirectionalFlow()
121+
inlet_face = ([open_boundary_size[1], 0.0], [open_boundary_size[1], channel.fluid_size[2]])
122+
inlet_reference_velocity = nothing
123+
inlet_reference_pressure = 0.2
124+
inlet_boundary_zone = BoundaryZone(; boundary_face=inlet_face, face_normal=flow_direction,
125+
open_boundary_layers, density=fluid_density,
126+
particle_spacing,
127+
reference_velocity=inlet_reference_velocity,
128+
reference_pressure=inlet_reference_pressure,
129+
initial_condition=inlet.fluid,
130+
boundary_type=inlet_boundary_type)
131+
132+
outlet_boundary_type = BidirectionalFlow()
133+
outlet_face = ([channel.fluid_size[1] - open_boundary_size[1], 0.0],
134+
[channel.fluid_size[1] - open_boundary_size[1], channel.fluid_size[2]])
135+
outlet_reference_velocity = nothing
136+
outlet_reference_pressure = 0.1
137+
outlet_flow_direction = (.-(flow_direction))
138+
outlet_boundary_zone = BoundaryZone(; boundary_face=outlet_face,
139+
face_normal=outlet_flow_direction,
140+
open_boundary_layers, density=fluid_density,
141+
particle_spacing,
142+
reference_velocity=outlet_reference_velocity,
143+
reference_pressure=outlet_reference_pressure,
144+
initial_condition=outlet.fluid,
145+
boundary_type=outlet_boundary_type)
146+
147+
open_boundary = OpenBoundarySystem(inlet_boundary_zone, outlet_boundary_zone; fluid_system,
138148
boundary_model=open_boundary_model,
139149
calculate_flow_rate=true,
140150
buffer_size=n_buffer_particles)
141151

142152
# ==========================================================================================
143153
# ==== Boundary
144-
wall = union(pipe.boundary)
154+
wall_boundary = union(channel.boundary)
145155

146-
boundary_model = BoundaryModelDummyParticles(wall.density, wall.mass,
156+
boundary_model = BoundaryModelDummyParticles(wall_boundary.density, wall_boundary.mass,
147157
AdamiPressureExtrapolation(),
148158
state_equation=state_equation,
149159
viscosity=viscosity,
150160
smoothing_kernel, smoothing_length)
151161

152-
boundary_system = WallBoundarySystem(wall, boundary_model)
162+
boundary_system = WallBoundarySystem(wall_boundary, boundary_model)
153163

154164
# ==========================================================================================
155165
# ==== Simulation
156-
min_corner = minimum(wall.coordinates .- 2 * particle_spacing, dims=2)
157-
max_corner = maximum(wall.coordinates .+ 2 * particle_spacing, dims=2)
166+
min_corner = minimum(wall_boundary.coordinates .- 2 * particle_spacing, dims=2)
167+
max_corner = maximum(wall_boundary.coordinates .+ 2 * particle_spacing, dims=2)
158168

159169
nhs = GridNeighborhoodSearch{2}(; cell_list=FullGridCellList(; min_corner, max_corner),
160170
update_strategy=ParallelUpdate())

‎examples/fluid/poiseuille_flow_3d.jl‎

Lines changed: 205 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,205 @@
1+
# ==========================================================================================
2+
# 3D Hagen-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 3D Hagen-Poiseuille flow simulation in a circular pipe
10+
# including open boundary conditions.
11+
# ==========================================================================================
12+
using TrixiParticles
13+
using OrdinaryDiffEq
14+
15+
# ==========================================================================================
16+
# ==== Resolution
17+
channel_length = 0.004
18+
channel_radius = 0.0005
19+
channel_diameter = 2 * channel_radius
20+
21+
particle_spacing_factor = 15
22+
particle_spacing = channel_diameter / particle_spacing_factor
23+
24+
# Make sure that the kernel support of fluid particles at a boundary is always fully sampled
25+
boundary_layers = 4
26+
27+
# Make sure that the kernel support of fluid particles at an open boundary is always
28+
# fully sampled.
29+
open_boundary_layers = 10
30+
31+
# ==========================================================================================
32+
# ==== Experiment Setup
33+
tspan = (0.0, 1.0)
34+
35+
fluid_density = 1000.0
36+
reynolds_number = 50
37+
imposed_pressure_drop = 0.1
38+
dynamic_viscosity = sqrt(fluid_density * channel_radius^2 * imposed_pressure_drop /
39+
(4 * reynolds_number))
40+
41+
v_max = channel_radius^2 * imposed_pressure_drop / (4 * dynamic_viscosity * channel_length)
42+
43+
# To accurately capture the initial transient this value has to be increased to around 100 which doubles the runtime
44+
sound_speed_factor = 50
45+
sound_speed = sound_speed_factor * v_max
46+
47+
flow_direction = (1.0, 0.0, 0.0)
48+
outlet_flow_direction = (-flow_direction[1], -flow_direction[2], -flow_direction[3])
49+
50+
n_particles_length = ceil(Int, channel_length / particle_spacing)
51+
52+
wall_cross_section = SphereShape(particle_spacing, channel_radius, (0.0, 0.0),
53+
fluid_density, n_layers=boundary_layers,
54+
layer_outwards=true, sphere_type=RoundSphere())
55+
56+
# Extend 2D coordinates to 3D by adding x-coordinates
57+
wall_cross_section_coordinates = hcat(particle_spacing / 2 *
58+
ones(nparticles(wall_cross_section)),
59+
wall_cross_section.coordinates')'
60+
61+
wall_boundary = extrude_geometry(wall_cross_section_coordinates; particle_spacing,
62+
direction=collect(flow_direction),
63+
density=fluid_density,
64+
n_extrude=n_particles_length)
65+
66+
fluid_cross_section = SphereShape(particle_spacing, channel_radius, (0.0, 0.0),
67+
fluid_density, sphere_type=RoundSphere())
68+
69+
# Extend 2D coordinates to 3D by adding x-coordinates
70+
fluid_cross_section_coordinates = hcat(particle_spacing / 2 *
71+
ones(nparticles(fluid_cross_section)),
72+
fluid_cross_section.coordinates')'
73+
74+
# The analytical solution uses the full inflow-to-outflow extent x in [0, channel_length].
75+
# Since the pressure is prescribed at the free surface of each `BoundaryZone` instead of at
76+
# its `boundary_face`, the inlet and outlet free surfaces lie at x = 0 and
77+
# x = channel_length (= L), respectively. Thus, the boundary zones are effectively part of
78+
# the computational domain in this setup.
79+
inlet_particles = extrude_geometry(fluid_cross_section_coordinates; particle_spacing,
80+
direction=collect(flow_direction),
81+
density=fluid_density,
82+
n_extrude=open_boundary_layers)
83+
84+
fluid_offset = [open_boundary_layers * particle_spacing, 0, 0]
85+
fluid_particles = extrude_geometry(fluid_cross_section_coordinates .+ fluid_offset;
86+
particle_spacing,
87+
direction=collect(flow_direction),
88+
density=fluid_density,
89+
n_extrude=(n_particles_length - 2 * open_boundary_layers))
90+
91+
outlet_offset = [(n_particles_length - open_boundary_layers) * particle_spacing, 0, 0]
92+
outlet_particles = extrude_geometry(fluid_cross_section_coordinates .+ outlet_offset;
93+
particle_spacing,
94+
direction=collect(flow_direction),
95+
n_extrude=open_boundary_layers,
96+
density=fluid_density)
97+
98+
n_buffer_particles = 10 * nparticles(fluid_cross_section)
99+
100+
# ==========================================================================================
101+
# ==== Fluid
102+
smoothing_length = 2 * particle_spacing
103+
smoothing_kernel = WendlandC2Kernel{3}()
104+
105+
fluid_density_calculator = ContinuityDensity()
106+
107+
kinematic_viscosity = dynamic_viscosity / fluid_density
108+
viscosity = ViscosityAdami(nu=kinematic_viscosity)
109+
110+
background_pressure = 7 * sound_speed_factor / 10 * fluid_density * v_max^2
111+
shifting_technique = TransportVelocityAdami(; background_pressure)
112+
113+
state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density,
114+
exponent=1)
115+
116+
fluid_system = WeaklyCompressibleSPHSystem(fluid_particles, fluid_density_calculator,
117+
state_equation, smoothing_kernel,
118+
buffer_size=n_buffer_particles,
119+
shifting_technique=shifting_technique,
120+
density_diffusion=DensityDiffusionMolteniColagrossi(delta=0.1),
121+
smoothing_length, viscosity=viscosity)
122+
123+
# ==========================================================================================
124+
# ==== Open Boundary
125+
open_boundary_model = BoundaryModelDynamicalPressureZhang()
126+
127+
inlet_boundary_type = BidirectionalFlow()
128+
inlet_face = ([
129+
open_boundary_layers * particle_spacing,
130+
-channel_diameter,
131+
-channel_diameter
132+
],
133+
[
134+
open_boundary_layers * particle_spacing,
135+
-channel_diameter,
136+
channel_diameter
137+
],
138+
[open_boundary_layers * particle_spacing, channel_diameter, channel_diameter])
139+
inlet_reference_velocity = nothing
140+
inlet_reference_pressure = 0.2
141+
142+
inlet_zone = BoundaryZone(; boundary_face=inlet_face, face_normal=flow_direction,
143+
open_boundary_layers, density=fluid_density, particle_spacing,
144+
reference_velocity=inlet_reference_velocity,
145+
reference_pressure=inlet_reference_pressure,
146+
initial_condition=inlet_particles,
147+
boundary_type=inlet_boundary_type)
148+
149+
outlet_boundary_type = BidirectionalFlow()
150+
outlet_face = ([outlet_offset[1], -channel_diameter, -channel_diameter],
151+
[outlet_offset[1], -channel_diameter, channel_diameter],
152+
[outlet_offset[1], channel_diameter, channel_diameter])
153+
outlet_reference_velocity = nothing
154+
outlet_reference_pressure = 0.1
155+
156+
outlet_zone = BoundaryZone(; boundary_face=outlet_face,
157+
face_normal=outlet_flow_direction,
158+
open_boundary_layers, density=fluid_density, particle_spacing,
159+
reference_velocity=outlet_reference_velocity,
160+
reference_pressure=outlet_reference_pressure,
161+
initial_condition=outlet_particles,
162+
boundary_type=outlet_boundary_type)
163+
164+
open_boundary = OpenBoundarySystem(inlet_zone, outlet_zone; fluid_system,
165+
boundary_model=open_boundary_model,
166+
buffer_size=n_buffer_particles)
167+
168+
# ==========================================================================================
169+
# ==== Boundary
170+
boundary_model = BoundaryModelDummyParticles(wall_boundary.density, wall_boundary.mass,
171+
AdamiPressureExtrapolation(),
172+
state_equation=state_equation,
173+
viscosity=viscosity,
174+
smoothing_kernel, smoothing_length)
175+
176+
boundary_system = WallBoundarySystem(wall_boundary, boundary_model)
177+
178+
# ==========================================================================================
179+
# ==== Simulation
180+
min_corner = minimum(wall_boundary.coordinates .- 2 * particle_spacing, dims=2)
181+
max_corner = maximum(wall_boundary.coordinates .+ 2 * particle_spacing, dims=2)
182+
183+
neighborhood_search = GridNeighborhoodSearch{3}(;
184+
cell_list=FullGridCellList(; min_corner,
185+
max_corner),
186+
update_strategy=ParallelUpdate())
187+
188+
semi = Semidiscretization(fluid_system, open_boundary,
189+
boundary_system,
190+
neighborhood_search=neighborhood_search,
191+
parallelization_backend=PolyesterBackend())
192+
193+
ode_problem = semidiscretize(semi, tspan)
194+
195+
info_callback = InfoCallback(interval=20)
196+
saving_callback = SolutionSavingCallback(dt=0.01, prefix="", output_directory="out")
197+
extra_callback = nothing
198+
199+
callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback(), extra_callback)
200+
201+
sol = solve(ode_problem, RDPK3SpFSAL35(),
202+
abstol=1e-6, # Default abstol is 1e-6 (may need tuning to prevent boundary penetration)
203+
reltol=1e-4, # Default reltol is 1e-3 (may need tuning to prevent boundary penetration)
204+
dtmax=1e-2, # Limit stepsize to prevent crashing
205+
save_everystep=false, callback=callbacks)

0 commit comments

Comments
 (0)