TrixiParticles.jlOpen/examples/fluid/dam_break_oil_film_2d.jl at 6a83b83ff03887124a13fc5d842697d2fde55cd4 · 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
# 2D dam break simulation with an oil film on top

using TrixiParticles
using OrdinaryDiffEq

# Size parameters
H = 0.6
W = 2 * H

gravity = 9.81
tspan = (0.0, 5.7 / sqrt(gravity))

# Resolution
fluid_particle_spacing = H / 60

# Numerical settings
smoothing_length = 3.5 * fluid_particle_spacing
sound_speed = 20 * sqrt(gravity * H)

# Physical values
nu_water = 8.9e-7
nu_oil = 6e-5
nu_ratio = nu_water / nu_oil

nu_sim_oil = max(0.01 * smoothing_length * sound_speed, nu_oil)
nu_sim_water = nu_ratio * nu_sim_oil

oil_viscosity = ViscosityMorris(nu=nu_sim_oil)

# TODO: broken if both systems use surface tension
trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "dam_break_2d.jl"),
              sol=nothing, fluid_particle_spacing=fluid_particle_spacing, tspan=tspan,
              viscosity_fluid=ViscosityMorris(nu=nu_sim_water),
              smoothing_length=smoothing_length,
              gravity=gravity, density_diffusion=nothing, sound_speed=sound_speed,
              prefix="", reference_particle_spacing=fluid_particle_spacing)

# ==========================================================================================
# ==== Setup oil layer

oil_size = (W, 0.1 * H)
oil_density = 700.0
oil_eos = StateEquationCole(; sound_speed, reference_density=oil_density, exponent=1,
                            clip_negative_pressure=false)

oil = RectangularShape(fluid_particle_spacing,
                       round.(Int, oil_size ./ fluid_particle_spacing),
                       zeros(length(oil_size)), density=oil_density)

# move on top of the water
for i in axes(oil.coordinates, 2)
    oil.coordinates[:, i] .+= [0.0, H]
end

oil_state_equation = StateEquationCole(; sound_speed, reference_density=oil_density,
                                       exponent=1, clip_negative_pressure=false)
oil_system = WeaklyCompressibleSPHSystem(oil, fluid_density_calculator,
                                         oil_eos, smoothing_kernel,
                                         smoothing_length, viscosity=oil_viscosity,
                                         acceleration=(0.0, -gravity),
                                         surface_tension=SurfaceTensionAkinci(surface_tension_coefficient=0.01),
                                         correction=AkinciFreeSurfaceCorrection(oil_density),
                                         reference_particle_spacing=fluid_particle_spacing)

# oil_system = WeaklyCompressibleSPHSystem(oil, fluid_density_calculator,
#                                          oil_eos, smoothing_kernel,
#                                          smoothing_length, viscosity=oil_viscosity,
#                                          acceleration=(0.0, -gravity),
#                                          surface_tension=SurfaceTensionMorris(surface_tension_coefficient=0.03),
#                                          reference_particle_spacing=fluid_particle_spacing)

# ==========================================================================================
# ==== Simulation
semi = Semidiscretization(fluid_system, oil_system, boundary_system,
                          neighborhood_search=GridNeighborhoodSearch{2}(update_strategy=nothing),
                          parallelization_backend=PolyesterBackend())
ode = semidiscretize(semi, tspan)

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # This is overwritten by the stepsize callback
            save_everystep=false, callback=callbacks);