|
1 | | -# 2D dam break simulation based on |
2 | | -# |
3 | | -# S. Marrone, M. Antuono, A. Colagrossi, G. Colicchio, D. le Touzé, G. Graziani. |
4 | | -# "δ-SPH model for simulating violent impact flows". |
5 | | -# In: Computer Methods in Applied Mechanics and Engineering, Volume 200, Issues 13–16 (2011), pages 1526–1542. |
6 | | -# https://doi.org/10.1016/J.CMA.2010.12.016 |
7 | | - |
| 1 | +# 2D dam break simulation using an ImplicitIncompressible SPH system |
8 | 2 | using TrixiParticles |
9 | | -using OrdinaryDiffEq |
10 | | - |
11 | | -# Size parameters |
12 | | -H = 0.8 |
13 | | -W = 2 * H |
14 | | - |
15 | | -# ========================================================================================== |
16 | | -# ==== Resolution |
17 | | -fluid_particle_spacing = H / 40 |
18 | | - |
19 | | -# Change spacing ratio to 3 and boundary layers to 1 when using Monaghan-Kajtar boundary model |
20 | | -boundary_layers = 4 |
21 | | -spacing_ratio = 1 |
22 | | - |
23 | | -boundary_particle_spacing = fluid_particle_spacing / spacing_ratio |
24 | | - |
25 | | -# ========================================================================================== |
26 | | -# ==== Experiment Setup |
27 | | -gravity = 9.81 |
28 | 3 |
|
29 | | -tspan = (0.0, 5.7 / sqrt(gravity)) |
| 4 | +# Load setup from dam break example |
| 5 | +trixi_include(@__MODULE__, |
| 6 | + joinpath(examples_dir(), "fluid", "dam_break_2d.jl"), |
| 7 | + sol=nothing, ode=nothing) |
30 | 8 |
|
31 | | -# Boundary geometry and initial fluid particle positions |
32 | | -initial_fluid_size = (W, H) |
33 | | -tank_size = (floor(5.366 * H / boundary_particle_spacing) * boundary_particle_spacing, 3.0) |
| 9 | +# Change smoohing kernel and length |
| 10 | +smoothing_length = 1.25 * fluid_particle_spacing |
| 11 | +smoothing_kernel = GaussianKernel{2}() |
34 | 12 |
|
35 | | -fluid_density = 1000.0 |
36 | | -sound_speed = 20 * sqrt(gravity * H) |
37 | | - |
38 | | -tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fluid_density, |
39 | | - n_layers=boundary_layers, spacing_ratio=spacing_ratio) |
40 | | - #acceleration=(0.0, -gravity)) |
41 | | - |
42 | | -# ========================================================================================== |
43 | | -# ==== Fluid |
44 | | -smoothing_length = 2.1 * fluid_particle_spacing |
45 | | -smoothing_kernel = WendlandC6Kernel{2}() |
46 | | - |
47 | | -# viscosity = ArtificialViscosityMonaghan(alpha=0.02, beta=0.0) |
| 13 | +# Calculate nu for the viscosity model |
48 | 14 | nu = 0.02 * smoothing_length * sound_speed/8 |
49 | | -# viscosity = ViscosityMorris(nu=nu) |
50 | | -viscosity = ViscosityAdami(nu=nu) |
51 | | -# Alternatively the density diffusion model by Molteni & Colagrossi can be used, |
52 | | -# which will run faster. |
53 | | -# density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1) |
54 | | -# density_diffusion = DensityDiffusionAntuono(tank.fluid, delta=0.1) |
55 | | -time_step=0.001 |
56 | | -fluid_system = ImplicitIncompressibleSPHSystem(tank.fluid, smoothing_kernel, |
57 | | - smoothing_length, fluid_density, |
58 | | - viscosity=viscosity, |
59 | | - acceleration=(0.0, -gravity), omega=0.5, |
60 | | - min_iterations=5, max_iterations=30, |
61 | | - time_step=time_step) |
62 | | - |
63 | | -# ========================================================================================== |
64 | | -# ==== Boundary |
65 | | -boundary_density_calculator = PressureMirroring() |
66 | | -boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, |
67 | | - state_equation=nothing, |
68 | | - boundary_density_calculator, |
69 | | - smoothing_kernel, smoothing_length) |
70 | 15 |
|
71 | | -boundary_system = BoundarySPHSystem(tank.boundary, boundary_model) |
72 | | - |
73 | | -# ========================================================================================== |
74 | | -# ==== Simulation |
75 | | -# `nothing` will automatically choose the best update strategy. This is only to be able |
76 | | -# to change this with `trixi_include`. |
77 | | -# semi = Semidiscretization(fluid_system, boundary_system, |
78 | | -# neighborhood_search=GridNeighborhoodSearch{2}(update_strategy=nothing), |
79 | | -# parallelization_backend=PolyesterBackend()) |
80 | | -semi = Semidiscretization(fluid_system, boundary_system) |
81 | | -ode = semidiscretize(semi, tspan) |
82 | | -#= |
83 | | -info_callback = InfoCallback(interval=100) |
84 | | -
|
85 | | -solution_prefix = "" |
86 | | -saving_callback = SolutionSavingCallback(dt=0.02, prefix=solution_prefix) |
87 | | -
|
88 | | -# Save at certain timepoints which allows comparison to the results of Marrone et al., |
89 | | -# i.e. (1.5, 2.36, 3.0, 5.7, 6.45). |
90 | | -# Please note that the images in Marrone et al. are obtained at a particle_spacing = H/320, |
91 | | -# which takes between 2 and 4 hours. |
92 | | -saving_paper = SolutionSavingCallback(save_times=[0.0, 0.371, 0.584, 0.743, 1.411, 1.597], |
93 | | - prefix="marrone_times") |
94 | | -
|
95 | | -# This can be overwritten with `trixi_include` |
96 | | -extra_callback = nothing |
97 | | -
|
98 | | -use_reinit = false |
99 | | -density_reinit_cb = use_reinit ? |
100 | | - DensityReinitializationCallback(semi.systems[1], interval=10) : |
101 | | - nothing |
102 | | -stepsize_callback = StepsizeCallback(cfl=0.9) |
103 | | -
|
104 | | -callbacks = CallbackSet(info_callback, saving_callback, stepsize_callback, extra_callback, |
105 | | - density_reinit_cb, saving_paper) |
106 | | -
|
107 | | -sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), |
108 | | - dt=1.0, # This is overwritten by the stepsize callback |
109 | | - save_everystep=false, callback=callbacks); |
110 | | -=# |
111 | | -info_callback = InfoCallback(interval=100) |
112 | | -saving_callback = SolutionSavingCallback(dt=0.02, prefix="") |
113 | | - |
114 | | -callbacks = CallbackSet(info_callback, saving_callback) |
115 | | - |
116 | | -sol = solve(ode, SymplecticEuler(), |
117 | | - dt=time_step, |
118 | | - save_everystep=false, callback=callbacks); |
119 | | - |
120 | | -plot(sol) |
| 16 | +# Use IISPH as fluid system |
| 17 | +time_step=0.001 |
| 18 | +min_iterations=10 |
| 19 | +max_iterations=30 |
| 20 | +IISPH_system = ImplicitIncompressibleSPHSystem(tank.fluid, smoothing_kernel, |
| 21 | + smoothing_length, fluid_density, |
| 22 | + viscosity=ViscosityAdami(nu=nu), |
| 23 | + acceleration=(0.0, -gravity), |
| 24 | + min_iterations=min_iterations, |
| 25 | + max_iterations=max_iterations, |
| 26 | + time_step=time_step) |
| 27 | + |
| 28 | +# Run the dam break simulation with these changes |
| 29 | +trixi_include(@__MODULE__, |
| 30 | + joinpath(examples_dir(), "fluid", "dam_break_2d.jl"), |
| 31 | + viscosity=ViscosityAdami(nu=nu), |
| 32 | + fluid_system=IISPH_system, |
| 33 | + boundary_density_calculator=PressureZeroing(), |
| 34 | + state_equation=nothing, |
| 35 | + callbacks=CallbackSet(info_callback, saving_callback), |
| 36 | + solver_function=SymplecticEuler(), dt=time_step) |
0 commit comments