add BoundaryModelTafuni

examples/fluid/pipe_flow_2d.jl

@@ -5,16 +5,16 @@ using OrdinaryDiffEq
 
 # ==========================================================================================
 # ==== Resolution
-particle_spacing = 0.05
+particle_spacing = 0.02
 
 # Make sure that the kernel support of fluid particles at a boundary is always fully sampled
-boundary_layers = 3
+boundary_layers = 5
 
 # Make sure that the kernel support of fluid particles at an open boundary is always
 # fully sampled.
 # Note: Due to the dynamics at the inlets and outlets of open boundaries,
 # it is recommended to use `open_boundary_layers > boundary_layers`
-open_boundary_layers = 6
+open_boundary_layers = 8
 
 # ==========================================================================================
 # ==== Experiment Setup
@@ -34,12 +34,11 @@ fluid_density = 1000.0
 
 # For this particular example, it is necessary to have a background pressure.
 # Otherwise the suction at the outflow is to big and the simulation becomes unstable.
-pressure = 1000.0
+pressure = 0.0
 
-sound_speed = 10 * prescribed_velocity
+sound_speed = 20 * prescribed_velocity
 
-state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density,
-                                   exponent=7, background_pressure=pressure)
+state_equation = nothing
 
 pipe = RectangularTank(particle_spacing, domain_size, boundary_size, fluid_density,
                        pressure=pressure, n_layers=boundary_layers,
@@ -48,7 +47,7 @@ pipe = RectangularTank(particle_spacing, domain_size, boundary_size, fluid_densi
 # Shift pipe walls in negative x-direction for the inflow
 pipe.boundary.coordinates[1, :] .-= particle_spacing * open_boundary_layers
 
-n_buffer_particles = 4 * pipe.n_particles_per_dimension[2]
+n_buffer_particles = 40 * pipe.n_particles_per_dimension[2]
 
 # ==========================================================================================
 # ==== Fluid
@@ -67,6 +66,8 @@ fluid_system = EntropicallyDampedSPHSystem(pipe.fluid, smoothing_kernel, smoothi
                                            buffer_size=n_buffer_particles)
 
 # Alternatively the WCSPH scheme can be used
+# state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density,
+#                                    exponent=1, background_pressure=pressure)
 # alpha = 8 * kinematic_viscosity / (smoothing_length * sound_speed)
 # viscosity = ArtificialViscosityMonaghan(; alpha, beta=0.0)
 
@@ -88,30 +89,27 @@ inflow = InFlow(; plane=([0.0, 0.0], [0.0, domain_size[2]]), flow_direction,
                 open_boundary_layers, density=fluid_density, particle_spacing)
 
 open_boundary_in = OpenBoundarySPHSystem(inflow; fluid_system,
-                                         boundary_model=BoundaryModelLastiwka(),
+                                         boundary_model=BoundaryModelTafuni(),
                                          buffer_size=n_buffer_particles,
-                                         reference_density=fluid_density,
-                                         reference_pressure=pressure,
                                          reference_velocity=velocity_function)
 
 outflow = OutFlow(; plane=([domain_size[1], 0.0], [domain_size[1], domain_size[2]]),
                   flow_direction, open_boundary_layers, density=fluid_density,
                   particle_spacing)
 
 open_boundary_out = OpenBoundarySPHSystem(outflow; fluid_system,
-                                          boundary_model=BoundaryModelLastiwka(),
+                                          boundary_model=BoundaryModelTafuni(),
                                           buffer_size=n_buffer_particles,
-                                          reference_density=fluid_density,
-                                          reference_pressure=pressure,
-                                          reference_velocity=velocity_function)
+                                          # reference_velocity=velocity_function,
+                                          reference_pressure=pressure)
 
 # ==========================================================================================
 # ==== Boundary
 
 boundary_model = BoundaryModelDummyParticles(pipe.boundary.density, pipe.boundary.mass,
                                              AdamiPressureExtrapolation(),
                                              state_equation=state_equation,
-                                             #viscosity=ViscosityAdami(nu=1e-4),
+                                             viscosity=ViscosityAdami(nu=1e-4),
                                              smoothing_kernel, smoothing_length)
 
 boundary_system = BoundarySPHSystem(pipe.boundary, boundary_model)

src/TrixiParticles.jl

@@ -69,7 +69,7 @@ export ArtificialViscosityMonaghan, ViscosityAdami, ViscosityMorris
 export DensityDiffusion, DensityDiffusionMolteniColagrossi, DensityDiffusionFerrari,
        DensityDiffusionAntuono
 export BoundaryModelMonaghanKajtar, BoundaryModelDummyParticles, AdamiPressureExtrapolation,
-       PressureMirroring, PressureZeroing, BoundaryModelLastiwka
+       PressureMirroring, PressureZeroing, BoundaryModelLastiwka, BoundaryModelTafuni
 export BoundaryMovement
 export examples_dir, validation_dir, trixi_include
 export trixi2vtk

src/schemes/boundary/boundary.jl

@@ -1,6 +1,7 @@
 include("dummy_particles/dummy_particles.jl")
 include("system.jl")
 include("open_boundary/boundary_zones.jl")
+include("open_boundary/mirroring.jl")
 include("open_boundary/method_of_characteristics.jl")
 include("open_boundary/system.jl")
 # Monaghan-Kajtar repulsive boundary particles require the `BoundarySPHSystem`

src/schemes/boundary/open_boundary/method_of_characteristics.jl

@@ -6,7 +6,12 @@ This model uses the characteristic variables to propagate the appropriate values
 to the outlet or inlet and have been proposed by Lastiwka et al. (2009). For more information
 about the method see [description below](@ref method_of_characteristics).
 """
-struct BoundaryModelLastiwka end
+struct BoundaryModelLastiwka
+    extrapolate_reference_values::Bool
+    function BoundaryModelLastiwka(; extrapolate_reference_values::Bool=true)
+        return new{}(extrapolate_reference_values)
+    end
+end
 
 @inline function update_quantities!(system, boundary_model::BoundaryModelLastiwka,
                                     v, u, v_ode, u_ode, semi, t)
@@ -15,6 +20,15 @@ struct BoundaryModelLastiwka end
 
     sound_speed = system_sound_speed(system.fluid_system)
 
+    if boundary_model.extrapolate_reference_values
+        (; prescribed_pressure, prescribed_velocity, prescribed_density) = cache
+
+        @trixi_timeit timer() "extrapolate and correct values" begin
+            extrapolate_values!(system, v_ode, u_ode, semi, t; prescribed_pressure,
+                                prescribed_velocity, prescribed_density)
+        end
+    end
+
     # Update quantities based on the characteristic variables
     @threaded system for particle in each_moving_particle(system)
         particle_position = current_coords(u, system, particle)
@@ -23,16 +37,16 @@ struct BoundaryModelLastiwka end
         J2 = cache.characteristics[2, particle]
         J3 = cache.characteristics[3, particle]
 
-        rho_ref = reference_value(reference_density, density[particle], system, particle,
+        rho_ref = reference_value(reference_density, density[particle],
                                   particle_position, t)
         density[particle] = rho_ref + ((-J1 + 0.5 * (J2 + J3)) / sound_speed^2)
 
-        p_ref = reference_value(reference_pressure, pressure[particle], system, particle,
+        p_ref = reference_value(reference_pressure, pressure[particle],
                                 particle_position, t)
         pressure[particle] = p_ref + 0.5 * (J2 + J3)
 
         v_current = current_velocity(v, system, particle)
-        v_ref = reference_value(reference_velocity, v_current, system, particle,
+        v_ref = reference_value(reference_velocity, v_current,
                                 particle_position, t)
         rho = density[particle]
         v_ = v_ref + ((J2 - J3) / (2 * sound_speed * rho)) * flow_direction
@@ -134,15 +148,16 @@ function evaluate_characteristics!(system, neighbor_system::FluidSystem,
 
         # Determine current and prescribed quantities
         rho_b = particle_density(v_neighbor_system, neighbor_system, neighbor)
-        rho_ref = reference_value(reference_density, density, system, particle,
+        rho_ref = reference_value(reference_density, density[particle],
                                   neighbor_position, t)
 
         p_b = particle_pressure(v_neighbor_system, neighbor_system, neighbor)
-        p_ref = reference_value(reference_pressure, pressure, system, particle,
+        p_ref = reference_value(reference_pressure, pressure[particle],
                                 neighbor_position, t)
 
         v_b = current_velocity(v_neighbor_system, neighbor_system, neighbor)
-        v_neighbor_ref = reference_value(reference_velocity, v, system, particle,
+        v_particle = current_velocity(v, system, particle)
+        v_neighbor_ref = reference_value(reference_velocity, v_particle,
                                          neighbor_position, t)
 
         # Determine characteristic variables

src/schemes/boundary/open_boundary/mirroring.jl

@@ -0,0 +1,182 @@
+# Two ways of providing the information to an open boundary are considered:
+# - physical quantities are either assigned a priori
+# - or extrapolated from the fluid domain to the buffer zones (inflow and outflow) using ghost nodes
+
+# The position of the ghost nodes is obtained by mirroring the boundary particles
+# into the fluid along a direction that is normal to the open boundary.
+
+# In order to calculate fluid quantities at the ghost nodes, a standard particle interpolation
+# would not be consistent due to the proximity of these points to an open boundary,
+# which translates into the lack of a full kernel support.
+struct BoundaryModelTafuni end
+
+function update_quantities!(system, ::BoundaryModelTafuni, v, u, v_ode, u_ode, semi, t)
+    @trixi_timeit timer() "extrapolate and correct values" begin
+        extrapolate_values!(system, v_ode, u_ode, semi, t; system.cache...)
+    end
+end
+
+update_final!(system, ::BoundaryModelTafuni, v, u, v_ode, u_ode, semi, t) = system
+
+function extrapolate_values!(system, v_ode, u_ode, semi, t; prescribed_density=false,
+                             prescribed_pressure=false, prescribed_velocity=false)
+    (; pressure, density, fluid_system, boundary_zone, reference_density,
+    reference_velocity, reference_pressure) = system
+
+    state_equation = system_state_equation(system.fluid_system)
+
+    # Static indices to avoid allocations
+    two_to_end = SVector{ndims(system)}(2:(ndims(system) + 1))
+
+    v_open_boundary = wrap_v(v_ode, system, semi)
+    v_fluid = wrap_v(v_ode, fluid_system, semi)
+    u_open_boundary = wrap_u(u_ode, system, semi)
+    u_fluid = wrap_u(u_ode, fluid_system, semi)
+
+    # Use the fluid-fluid nhs, since the boundary particles are mirrored into the fluid domain
+    neighborhood_search = get_neighborhood_search(fluid_system, fluid_system, semi)
+
+    @threaded system for particle in each_moving_particle(system)
+        particle_coords = current_coords(u_open_boundary, system, particle)
+        ghost_node_position = mirror_position(particle_coords, boundary_zone)
+
+        # Set zero
+        correction_matrix = zero(SMatrix{ndims(system) + 1, ndims(system) + 1,
+                                         eltype(system)})
+        interpolated_pressure = zero(SVector{ndims(system) + 1, eltype(system)})
+
+        interpolated_velocity = zero(SMatrix{ndims(system), ndims(system) + 1,
+                                             eltype(system)})
+
+        for neighbor in PointNeighbors.eachneighbor(ghost_node_position,
+                                                    neighborhood_search)
+            neighbor_coords = current_coords(u_fluid, fluid_system, neighbor)
+            pos_diff = ghost_node_position - neighbor_coords
+            distance2 = dot(pos_diff, pos_diff)
+
+            if distance2 <= neighborhood_search.search_radius^2
+                distance = sqrt(distance2)
+
+                m_b = hydrodynamic_mass(fluid_system, neighbor)
+                rho_b = particle_density(v_fluid, fluid_system, neighbor)
+                pressure_b = particle_pressure(v_fluid, fluid_system, neighbor)
+                v_b = current_velocity(v_fluid, fluid_system, neighbor)
+
+                # Project `v_b` to the normal direction of the boundary zone
+                # See https://doi.org/10.1016/j.jcp.2020.110029 Section 3.3.:
+                # "Because flow from the inlet interface occurs perpendicular to the boundary,
+                # only this component of interpolated velocity is kept [...]"
+                v_b = dot(v_b, boundary_zone.flow_direction) * boundary_zone.flow_direction
+
+                kernel_value = smoothing_kernel(fluid_system, distance)
+                grad_kernel = smoothing_kernel_grad(fluid_system, pos_diff, distance)
+
+                L, R = correction_arrays(kernel_value, grad_kernel, pos_diff, rho_b, m_b)
+
+                correction_matrix += L
+
+                !(prescribed_pressure) && (interpolated_pressure += pressure_b * R)
+
+                !(prescribed_velocity) && (interpolated_velocity += v_b * R')
+            end
+        end
+
+        if abs(det(correction_matrix)) < 1e-9
+            L_inv = I
+        else
+            L_inv = inv(correction_matrix)
+        end
+
+        pos_diff = particle_coords - ghost_node_position
+
+        # In Negi et al. (2020) https://doi.org/10.1016/j.cma.2020.113119,
+        # they have modified the equation for extrapolation to
+        #
+        #           f_0 = f_k - (r_0 - r_k) ⋅ ∇f_k
+        #
+        # in order to get zero gradient at the outlet interface.
+        if prescribed_pressure
+            pressure[particle] = reference_value(reference_pressure, pressure[particle],
+                                                 particle_coords, t)
+        else
+            f_p = L_inv * interpolated_pressure
+            df_p = f_p[two_to_end]
+
+            pressure[particle] = f_p[1] + dot(pos_diff, df_p)
+        end
+
+        if prescribed_velocity
+            v_particle = current_velocity(v_open_boundary, system, particle)
+            v_ref = reference_value(reference_velocity, v_particle, particle_coords, t)
+            @inbounds for dim in 1:ndims(system)
+                v_open_boundary[dim, particle] = v_ref[dim]
+            end
+        else
+            @inbounds for dim in 1:ndims(system)
+                f_v = L_inv * interpolated_velocity[dim, :]
+                df_v = f_v[two_to_end]
+
+                v_open_boundary[dim, particle] = f_v[1] + dot(pos_diff, df_v)
+            end
+        end
+
+        # Unlike Tafuni et al. (2018), we calculate the density using the inverse state-equation.
+        if prescribed_density
+            density[particle] = reference_value(reference_density, density[particle],
+                                                particle_coords, t)
+        else
+            inverse_state_equation!(density, state_equation, pressure, particle)
+        end
+    end
+
+    return system
+end
+
+function correction_arrays(W_ab, grad_W_ab, pos_diff::SVector{3}, rho_b, m_b)
+    x_ab = pos_diff[1]
+    y_ab = pos_diff[2]
+    z_ab = pos_diff[3]
+
+    grad_W_ab_x = grad_W_ab[1]
+    grad_W_ab_y = grad_W_ab[2]
+    grad_W_ab_z = grad_W_ab[3]
+
+    V_b = m_b / rho_b
+
+    M_ = @SMatrix [W_ab W_ab*x_ab W_ab*y_ab W_ab*z_ab;
+                   grad_W_ab_x grad_W_ab_x*x_ab grad_W_ab_x*y_ab grad_W_ab_x*z_ab;
+                   grad_W_ab_y grad_W_ab_y*x_ab grad_W_ab_y*y_ab grad_W_ab_y*z_ab;
+                   grad_W_ab_z grad_W_ab_z*x_ab grad_W_ab_z*y_ab grad_W_ab_z*z_ab]
+
+    L = V_b * M_
+
+    R = V_b * SVector(W_ab, grad_W_ab_x, grad_W_ab_y, grad_W_ab_z)
+
+    return L, R
+end
+
+function correction_arrays(W_ab, grad_W_ab, pos_diff::SVector{2}, rho_b, m_b)
+    x_ab = pos_diff[1]
+    y_ab = pos_diff[2]
+
+    grad_W_ab_x = grad_W_ab[1]
+    grad_W_ab_y = grad_W_ab[2]
+
+    V_b = m_b / rho_b
+
+    M_ = @SMatrix [W_ab W_ab*x_ab W_ab*y_ab;
+                   grad_W_ab_x grad_W_ab_x*x_ab grad_W_ab_x*y_ab;
+                   grad_W_ab_y grad_W_ab_y*x_ab grad_W_ab_y*y_ab]
+    L = V_b * M_
+
+    R = V_b * SVector(W_ab, grad_W_ab_x, grad_W_ab_y)
+
+    return L, R
+end
+
+function mirror_position(particle_coords, boundary_zone)
+    particle_position = particle_coords - boundary_zone.zone_origin
+    dist = dot(particle_position, boundary_zone.flow_direction)
+
+    return particle_coords - 2 * dist * boundary_zone.flow_direction
+end