Download this testcase.

Water–Salt Water Avalanche¶

This test case simulates the collapse of a dense granular column immersed in a two-fluid environment composed of fresh water and sea water. The grain column is initially surrounded by denser sea water, while the rest of the domain is filled with lighter fresh water.

The case cannot be run standalone: it requires the output of the depot testcase (../depot/depot.py) to initialize a compact granular column. The script can, however, generate a synthetic initial packing using a hexagonal arrangement for testing purposes.

Keywords¶

FEM, DEM, Unresolved, Friction, Two-Fluids

import sys, os, shutil, subprocess
from migflow import fluid, scontact, time_integration
import numpy as np
import time

Output Directory¶

Create an output directory for results and generate a 2D mesh using Gmsh.

outputdir = "output" if len(sys.argv) < 2 else sys.argv[1]
geomesh_filename = f"mesh_two_fluids.geo"
mesh_filename = f"{outputdir}/mesh.msh"
shutil.rmtree(outputdir, ignore_errors=True)
os.makedirs(outputdir)
subprocess.call(["gmsh", "-2", geomesh_filename, "-o", mesh_filename])

Particle Packing Function¶

Define a helper function to create an initial hexagonal particle arrangement when no pre-deposited configuration is available. The packing density is controlled by the range of particle radii.

use_pre_deposit = False


def hexagonal_packing(p, x0, y0, lx, ly, rmin, rmax):
    """Set all the particles’ center positions using a hexagonal arrangement.

    Parameters
    ----------
    p : migflow.scontact.ParticleProblem
        Particle problem object.
    x0, y0 : float
        Coordinates of the bottom-left corner of the packing region.
    lx, ly : float
        Width and height of the region to fill with particles.
    rmin, rmax : float
        Minimum and maximum particle radii.
    """
    x = np.arange(x0 + rmax, x0 + lx - rmax, 2 * rmax)
    y = np.arange(y0 + rmax, y0 + ly - rmax, 2 * (3**0.5) * rmax)
    x, y = np.meshgrid(x, y)
    x1 = np.arange(2 * rmax + x0, x0 + lx - rmax, 2 * rmax)
    y1 = np.arange(y0 + 3**0.5 * rmax + rmax, y0 + ly - rmax, 2 * (3**0.5) * rmax)
    x1, y1 = np.meshgrid(x1, y1)
    x = np.concatenate([x.reshape(-1), x1.reshape(-1)])
    y = np.concatenate([y.reshape(-1), y1.reshape(-1)])
    order = np.argsort(y)
    x = x[order]
    y = y[order]
    for xi, yi in zip(x, y):
        z = np.random.random((1))
        r = rmin + (rmax - rmin) * z
        p.add_particle((xi, yi), r, r**2 * np.pi * rhop, material="particle")

Physical Parameters¶

Define the physical properties for the particles, fresh water, and sea water. The density contrast between the two fluids drives buoyancy and stratification.

g = np.array([0, -9.81])  # gravity vector [m/s²]
rhop = 1500  # particle density [kg/m³]
r = 2e-3  # particle radius [m]
# Fresh water (light phase)
rhof = 1000  # density [kg/m³]
nuf = 1e-6  # kinematic viscosity [m²/s]
# Sea water (dense phase)
rhom = 1050  # density [kg/m³]
num = 1e-6  # kinematic viscosity [m²/s]

Particle Problem¶

Create the particle problem, load boundaries from the mesh, and initialize the grain assembly either from a pre-deposited configuration or by generating a regular packing.

p = scontact.ParticleProblem(2)
p.set_fixed_contact_geometry(0)  # workaround for contact geometry flag
p.load_msh_boundaries(
    mesh_filename, ["Top", "Bottom", "Left", "Right"], material="Wall"
)
if use_pre_deposit:
    p1 = scontact.ParticleProblem(2)
    p1.read_mig("../depot/output", iteration=-1)
    p.add_particles(p1.position(), p1.r(), p1.mass(), v=p1.velocity())
else:
    hexagonal_packing(p, 0, 0, 0.1, 0.199, 0.9 * r, 1.1 * r)
p.set_friction_coefficient(1.0, "Wall", "Particle")
p.set_friction_coefficient(0.5, "Particle", "Particle")

Fluid Problem¶

Define the fluid problem with two phases (fresh water and sea water). The density and viscosity fields are specified as lists for the two components. The particle information (positions, velocities, contact forces) is passed to the fluid solver for coupling.

f = fluid.FluidProblem(2, g, [nuf * rhof, num * rhom], [rhof, rhom])
f.load_msh(mesh_filename)
f.set_wall_boundary("Bottom")
f.set_wall_boundary("Left")
f.set_wall_boundary("Right")
f.set_wall_boundary("Top")
f.set_mean_pressure(0)
f.set_particles(
    p.delassus(), p.volume(), p.position(), p.velocity(), p.omega(), p.contact_forces()
)

Initial Concentration Field¶

The concentration field defines the spatial distribution of sea water (c = 1) and fresh water (c = 0). Initially, the left region is filled with sea water to represent the salt water column surrounding the grains.

x = f.coordinates()
c = np.zeros(f.n_nodes())
c[x[:, 0] < 0.1] = 1
f.set_concentration_cg(c)

Time Integration¶

The simulation runs for a short time interval with a fixed time step. At each iteration, the fluid phase is solved and the fluid-grain forces are applied to the particles, and then system evolves according to the contact dynamics solver. The simulation runs until the avalanche stabilizes.

t = 0  # initial time
i = 0  # iteration counter
tic = time.time()
outf = 10  # number of iterations between outputs
dt = 1e-3  # time step [s]
tEnd = 2  # final time [s]
mass = np.pi * p.r() ** 2 * rhop  # particle masses
while t < tEnd:
    print(f"{i = :4d} ----- {t = :g}")
    if i % outf == 0:
        p.write_mig(outputdir, t)
        f.write_mig(outputdir, t)
    time_integration.iterate(f, p, dt, min_nsub=10, external_particles_forces=g * mass)
    t += dt
    i += 1

Plot¶

python3 -m migflow.plot.migplot output --actors fluid particles --fluid-field concentration