TrixiParticles.jlOpen/examples/fsi/hydrostatic_water_column_2d.jl at 2966228d6c428afcad80c4d320244795dd00b3a9 · svchb/TrixiParticles.jlOpen · GitHub

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
147
148
149
150
151
152
153
154
155
156
157
158
159
160
# ==========================================================================================
# 2D Hydrostatic Water Column on an Elastic Plate (FSI)
#
# Case "Elastic plate under a hydrostatic water column" as described in
# "A fluid–structure interaction model for free-surface flows and flexible structures
# using smoothed particle hydrodynamics on a GPU" by J. O'Connor and B.D. Rogers
# published in Journal of Fluids and Structures
# https://doi.org/10.1016/j.jfluidstructs.2021.103312
# ==========================================================================================

using TrixiParticles
using OrdinaryDiffEq
using JSON

# ==========================================================================================
# ==== Options
use_edac = true  # Use EDAC or WCSPH

# ==========================================================================================
# ==== Experiment Parameters
# The paper uses 5 to 40 (3 for short runtime in CI)
n_particles_plate_y = 5
boundary_layers = 3
spacing_ratio = 1
gravity = 9.81
# Recommended: run at least to 0.5 (paper runs to 1.0); using 0.3 for speed in CI
tspan = (0.0, 0.5)

initial_fluid_size = (1.0, 2.0)
plate_size = (1.0, 0.05)

fluid_density = 1000.0
structure_density = 2700.0

E = 67.5e9
nu = 0.3
sound_speed = 50.0

# Particle spacings (structure and fluid share same spacing)
structure_particle_spacing = plate_size[2] / (n_particles_plate_y - 1)
fluid_particle_spacing = structure_particle_spacing
boundary_particle_spacing = fluid_particle_spacing / spacing_ratio

# Analytical solution (constant)
D = E * plate_size[2]^3 / (12 * (1 - nu^2))
analytical_value = -0.0026 * gravity *
                   (fluid_density * initial_fluid_size[2] +
                    structure_density * plate_size[2]) / D

# ==========================================================================================
# ==== Geometry Definitions: Tank, Beam, and Fixed Particles
n_particles_plate_x = round(Int, plate_size[1] / structure_particle_spacing + 1)
n_particles_per_dimension = (n_particles_plate_x, n_particles_plate_y)

plate = RectangularShape(structure_particle_spacing, n_particles_per_dimension,
                         (0.0, -plate_size[2]), density=structure_density,
                         place_on_shell=true)

left_wall = RectangularShape(structure_particle_spacing, (3, n_particles_plate_y),
                             (-3 * structure_particle_spacing, -plate_size[2]),
                             density=structure_density, place_on_shell=true)
right_wall = RectangularShape(structure_particle_spacing, (3, n_particles_plate_y),
                              (plate_size[1] + structure_particle_spacing,
                               -plate_size[2]), density=structure_density,
                              place_on_shell=true)
fixed_particles = union(left_wall, right_wall)

structure_geometry = union(plate, fixed_particles)

# ==========================================================================================
# ==== Smoothing Kernel, Boundary, and Related Quantities
smoothing_kernel = WendlandC2Kernel{2}()
smoothing_length_structure = sqrt(2) * structure_particle_spacing

smoothing_length_fluid = sqrt(2) * fluid_particle_spacing

hydrodynamic_densities = fluid_density * ones(size(structure_geometry.density))
hydrodynamic_masses = hydrodynamic_densities *
                      structure_particle_spacing^ndims(structure_geometry)
boundary_density_calculator = AdamiPressureExtrapolation()

# ==========================================================================================
# ==== Run Simulations
state_equation = use_edac ? nothing :
                 StateEquationCole(; sound_speed, reference_density=fluid_density,
                                   exponent=7, clip_negative_pressure=false)

tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, (plate_size[1], 3.0),
                       min_coordinates=(0.0, fluid_particle_spacing / 2),
                       fluid_density, n_layers=boundary_layers,
                       spacing_ratio=spacing_ratio,
                       faces=(true, true, false, false),
                       acceleration=(0.0, -gravity),
                       state_equation=state_equation)

if use_edac
    fluid_system = EntropicallyDampedSPHSystem(tank.fluid, smoothing_kernel,
                                               smoothing_length_fluid, sound_speed,
                                               acceleration=(0.0, -gravity),
                                               correction=ShepardKernelCorrection(),
                                               source_terms=SourceTermDamping(;
                                                                              damping_coefficient=0.05))
else
    fluid_density_calculator = ContinuityDensity()
    density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1)
    # density_diffusion = DensityDiffusionAntuono(tank.fluid, delta=0.1)
    fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator,
                                               state_equation, smoothing_kernel,
                                               smoothing_length_fluid,
                                               density_diffusion=density_diffusion,
                                               acceleration=(0.0, -gravity),
                                               source_terms=SourceTermDamping(;
                                                                              damping_coefficient=0.05))
end

boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass,
                                             state_equation=state_equation,
                                             boundary_density_calculator,
                                             smoothing_kernel, smoothing_length_fluid)
boundary_system = WallBoundarySystem(tank.boundary, boundary_model)
boundary_model_structure = BoundaryModelDummyParticles(hydrodynamic_densities,
                                                       hydrodynamic_masses,
                                                       state_equation=state_equation,
                                                       boundary_density_calculator,
                                                       smoothing_kernel,
                                                       smoothing_length_structure)
structure_system = TotalLagrangianSPHSystem(structure_geometry, smoothing_kernel,
                                            smoothing_length_structure,
                                            E, nu, boundary_model=boundary_model_structure,
                                            n_clamped_particles=nparticles(fixed_particles),
                                            acceleration=(0.0, -gravity))

min_corner = min.(minimum(structure_geometry.coordinates, dims=2),
                  minimum(tank.boundary.coordinates, dims=2),
                  minimum(tank.fluid.coordinates, dims=2)) .- smoothing_length_fluid
max_corner = max.(maximum(structure_geometry.coordinates, dims=2),
                  maximum(tank.boundary.coordinates, dims=2),
                  maximum(tank.fluid.coordinates, dims=2)) .+ smoothing_length_fluid

cell_list = FullGridCellList(; min_corner, max_corner)
neighborhood_search = GridNeighborhoodSearch{2}(; update_strategy=ParallelUpdate(),
                                                cell_list)
semi = Semidiscretization(structure_system, fluid_system, boundary_system,
                          neighborhood_search=neighborhood_search,
                          parallelization_backend=PolyesterBackend())
ode = semidiscretize(semi, tspan)

split_integration = SplitIntegrationCallback(CarpenterKennedy2N54(williamson_condition=false),
                                             callback=StepsizeCallback(cfl=1.6),
                                             dt=1.0) # dt will be overwritten by the stepsize cb

# This can be overwritten with `trixi_include`
extra_callback = nothing

info_callback = InfoCallback(interval=100)
saving_callback = SolutionSavingCallback(dt=0.1, prefix="")
callbacks = CallbackSet(info_callback, saving_callback, split_integration, extra_callback)

sol = solve(ode, RDPK3SpFSAL49(), dt=1e-8, reltol=1e-6, save_everystep=false,
            callback=callbacks)