-
Notifications
You must be signed in to change notification settings - Fork 21
Expand file tree
/
Copy pathfalling_spheres_2d.jl
More file actions
146 lines (119 loc) · 7.1 KB
/
Copy pathfalling_spheres_2d.jl
File metadata and controls
146 lines (119 loc) · 7.1 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
# ==========================================================================================
# 2D Falling Spheres in Fluid (FSI) - Base Setup
#
# This file provides a base setup for simulating one or two elastic spheres
# falling into a fluid in a tank.
# ==========================================================================================
using TrixiParticles
using OrdinaryDiffEqLowStorageRK
# ==========================================================================================
# ==== Resolution
fluid_particle_spacing = 0.02
structure_particle_spacing = fluid_particle_spacing
# Change spacing ratio to 3 and boundary layers to 1 when using Monaghan-Kajtar boundary model
boundary_layers = 3
spacing_ratio = 1
# ==========================================================================================
# ==== Experiment Setup
gravity = 9.81
tspan = (0.0, 1.0)
# Boundary geometry and initial fluid particle positions
initial_fluid_size = (2.0, 0.9)
tank_size = (2.0, 1.0)
fluid_density = 1000.0
sound_speed = 10 * sqrt(gravity * initial_fluid_size[2])
state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density,
exponent=1)
tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fluid_density;
n_layers=boundary_layers, spacing_ratio,
faces=(true, true, true, false), acceleration=(0.0, -gravity),
state_equation)
sphere1_radius = 0.3
sphere2_radius = 0.2
sphere1_density = 500.0
sphere2_density = 1100.0
# Young's modulus and Poisson ratio
sphere1_E = 7e4
sphere2_E = 1e5
nu = 0.0
sphere1_center = (0.5, 1.6)
sphere2_center = (1.5, 1.6)
sphere1 = SphereShape(structure_particle_spacing, sphere1_radius, sphere1_center,
sphere1_density, sphere_type=VoxelSphere())
sphere2 = SphereShape(structure_particle_spacing, sphere2_radius, sphere2_center,
sphere2_density, sphere_type=VoxelSphere())
# ==========================================================================================
# ==== Fluid
fluid_smoothing_length = 1.5 * fluid_particle_spacing
fluid_smoothing_kernel = WendlandC2Kernel{2}()
fluid_density_calculator = ContinuityDensity()
viscosity = ArtificialViscosityMonaghan(alpha=0.02, beta=0.0)
density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1)
fluid_system = WeaklyCompressibleSPHSystem(tank.fluid;
smoothing_kernel=fluid_smoothing_kernel,
smoothing_length=fluid_smoothing_length,
density_calculator=fluid_density_calculator,
state_equation, viscosity, density_diffusion,
acceleration=(0.0, -gravity))
# ==========================================================================================
# ==== Boundary
boundary_density_calculator = BernoulliPressureExtrapolation()
# Clip negative boundary pressure values to avoid sticking artifacts at the boundary.
boundary_model = BoundaryModelDummyParticles(tank.boundary; fluid_system=fluid_system,
boundary_density_calculator,
clip_negative_pressure=true)
boundary_system = WallBoundarySystem(tank.boundary, boundary_model)
# ==========================================================================================
# ==== Structure
structure_smoothing_length = sqrt(2) * structure_particle_spacing
structure_smoothing_kernel = WendlandC2Kernel{2}()
# For the FSI we need the hydrodynamic masses and densities in the structure boundary model
hydrodynamic_densites_1 = fluid_density * ones(size(sphere1.density))
hydrodynamic_masses_1 = hydrodynamic_densites_1 *
structure_particle_spacing^ndims(fluid_system)
structure_boundary_model_1 = BoundaryModelDummyParticles(hydrodynamic_densites_1,
hydrodynamic_masses_1,
boundary_density_calculator,
fluid_smoothing_kernel,
fluid_smoothing_length;
state_equation,
clip_negative_pressure=true)
hydrodynamic_densites_2 = fluid_density * ones(size(sphere2.density))
hydrodynamic_masses_2 = hydrodynamic_densites_2 *
structure_particle_spacing^ndims(fluid_system)
structure_boundary_model_2 = BoundaryModelDummyParticles(hydrodynamic_densites_2,
hydrodynamic_masses_2,
boundary_density_calculator,
fluid_smoothing_kernel,
fluid_smoothing_length;
state_equation,
clip_negative_pressure=true)
structure_system_1 = TotalLagrangianSPHSystem(sphere1;
smoothing_kernel=structure_smoothing_kernel,
smoothing_length=structure_smoothing_length,
young_modulus=sphere1_E,
poisson_ratio=nu,
acceleration=(0.0, -gravity),
boundary_model=structure_boundary_model_1,
penalty_force=PenaltyForceGanzenmueller(alpha=0.3))
structure_system_2 = TotalLagrangianSPHSystem(sphere2;
smoothing_kernel=structure_smoothing_kernel,
smoothing_length=structure_smoothing_length,
young_modulus=sphere2_E,
poisson_ratio=nu,
acceleration=(0.0, -gravity),
boundary_model=structure_boundary_model_2,
penalty_force=PenaltyForceGanzenmueller(alpha=0.3))
# ==========================================================================================
# ==== Simulation
semi = Semidiscretization(fluid_system, boundary_system, structure_system_1,
structure_system_2)
ode = semidiscretize(semi, tspan)
info_callback = InfoCallback(interval=50)
saving_callback = SolutionSavingCallback(dt=0.02, output_directory="out", prefix="")
callbacks = CallbackSet(info_callback, saving_callback)
# Use a Runge-Kutta method with automatic (error based) time step size control.
sol = solve(ode, RDPK3SpFSAL49(),
abstol=1e-6, # Default abstol is 1e-6
reltol=1e-3, # Default reltol is 1e-3
save_everystep=false, callback=callbacks);