# %% # Drag on a Fixed Disk in a 2D Immersed Granular Flow # =================================================== # This example studies the drag force acting on a fixed circular obstacle # (a “disk”) immersed in a granular flow saturated by a viscous fluid. # # %% # Keywords # -------- # Drag, DEM, FEM, Boundary Force, 2D-3D Ratio # # %% # Description # ----------- # The test demonstrates: # - How to set up a mixed granular–fluid simulation in a 2D geometry. # - How to impose an inflow boundary condition to simulate a stream. # - How to insert, remove, and constrain particles dynamically. # - How to compute and record the total drag force acting on a solid obstacle. # A 2D–3D scaling factor is applied to account for the reduced dimensionality. from migflow import scontact, fluid, time_integration, gmsh_io, volume_coupling as vc import numpy as np import os, sys, shutil, subprocess # %% # Output Directory and Mesh Generation # ------------------------------------ # The output directory is created, and a 2D annular mesh is generated using Gmsh. outputdir = "output_2d_drag" if len(sys.argv) < 2 else sys.argv[1] geomesh_filename = "2d_mesh.geo" mesh_filename = f"{outputdir}/mesh.msh" csv_file = f"{outputdir}/Drag.csv" shutil.rmtree(outputdir, ignore_errors=True) os.makedirs(outputdir) subprocess.call(["gmsh", "-2", geomesh_filename, "-o", mesh_filename]) # %% # Geometrical Parameters # ---------------------- h_p = 0.3 # particle bed height [m] w_p = 0.3 # particle bed width [m] w = 0.3 # particle bed width [m] h = 1.4 # total domain height [m] r_inner = 0.025 # radius of the fixed disk [m] ratio_2d_3d = 0.7 # 2D-3D correction # %% # Physical Parameters # ------------------- v_imposed = -0.05 # stream velocity rhop = 2500 # particle density rho = 1000 # fluid density r = 3e-3 # particle radius g = np.array([0, -9.81]) # gravity nu = 1e-6 # fluid kinematic viscosity mu = rho * nu friction = 0.2 # %% # Particle Problem Initialization # ------------------------------- # Load mesh boundaries and define contact properties. p = scontact.ParticleProblem(2) gmsh_io.load_msh_boundaries(p, mesh_filename, ["Inner", "Left", "Right"], material="Steel") p.set_friction_coefficient(friction, "Sand", "Sand") p.set_friction_coefficient(friction, "Sand", "Steel") # %% # Particle Generation # ------------------- # The following function generates particle coordinates on a regular grid # over a rectangular region. def gen_rect(origin, w, h, step): """Generate coordinates on a rectangular grid Keyword arguments: origin -- origin of top left corner w : width h : height step : maximal radius """ eps = 1e-4 * step x = np.arange(-w / 2 + step - eps, w / 2 - step + eps, 2 * step) + origin[0] y = np.arange(step, h - step, 2 * step) + origin[1] x, y = np.meshgrid(x, y) return x.reshape(-1), y.reshape(-1) x, y = gen_rect([0, -h / 2], w, 4 * h_p, r) for xi, yi in zip(x, y): if xi**2 + yi**2 > (r_inner + r) ** 2: p.add_particle((xi, yi), r, r**2 * np.pi * rhop, "Sand") x, y = gen_rect([0, 0.5 * h / 1.4], w, 0.5 * h_p, 1.2 * r) for xi, yi in zip(x, y): p.add_particle((xi, yi), r, r**2 * np.pi * rhop, "Sand") # %% # Fluid Problem Initialization # ---------------------------- # The fluid problem is loaded from the same mesh and coupled to the particle phase. f = fluid.FluidProblem2(g, nu * rho, rho) gmsh_io.load_msh(f, mesh_filename) f.set_wall_boundary("Inner") f.set_wall_boundary("Left", velocity=[0, v_imposed]) f.set_wall_boundary("Right", velocity=[0, v_imposed]) f.set_open_boundary("Bottom", velocity=[0, v_imposed]) f.set_strong_boundary("Left", velocity=[0, v_imposed]) f.set_strong_boundary("Right", velocity=[0, v_imposed]) f.set_strong_boundary("Bottom", velocity=[0, v_imposed]) f.set_open_boundary("Top", pressure=0) # %% # Numerical Parameters # -------------------- dt = 5e-3 # time step t = 0 # initial time tEnd = 1.5 # final time i = 0 # iteration number outf = 2 # iterations between data frames # %% # Computation Loop # ---------------- # Accumulate the CONTACT force the grains apply to the inner disk, averaged over # the DEM sub-steps. get_boundary_forces already returns the total for the named # boundary, shape (dim,) -- the old np.sum(..., axis=0) collapsed that to a # SCALAR, which then broadcast the same number into both components of F. It # never showed because the callback was not being invoked at all (iterate() was # dropping after_sub_iter on this branch), so both defects sat here unseen. def accumulate(F, n_divide): F += np.asarray(p.get_boundary_forces("Inner"), float).ravel() / n_divide while t < tEnd: print(f"{i:4d}, {t:.6g}/{tEnd:.6g}") # Remove particles that have fallen outside the domain alert = 4 * r keep = (p.body_position()[:, 1] > -h / 1.4 * 0.55) & ( p.body_position()[:, 1] < h / 2 - alert ) p.remove_bodies_flag(keep) # Constrain particles at the bottom nonzero = p.r()[:, 0] != 0 position = p.position()[nonzero, :] p.body_invert_mass()[p.body_position()[:, 1] < -0.5] = 0 p.body_invert_inertia()[p.body_position()[:, 1] < -0.5] = 0 a = p.body_velocity()[p.n_bodies() - position.shape[0] :, :] b = p.body_omega()[p.n_bodies() - position.shape[0] :, :] c = p.body_invert_mass()[p.n_bodies() - position.shape[0] :, :] a[(c == 0)[:, 0], :] = [0, v_imposed] b[(c == 0)[:, 0], 0] = 0 # Insert new particles at the top when needed if np.amax(position[:, 1]) + r < 0.5: for xi, yi in zip(x, y): p.add_particle((xi, yi), r, r**2 * np.pi * rhop, "Sand") # Set particles to the FEM module. The 2d-3d correction scales the body # volume, so tiers A and B see the corrected body: its radius and its # density both come from that volume. No datum -- this only refreshes the # coupling so the frame written below carries a porosity field, and # time_integration.iterate re-sets all three tiers, datum included. vol = p.volume() * ratio_2d_3d geo = vc.get_particles(f, p.position(), vc.body_radii(vol, 2), density=vc.body_density(p, volume=vol)) clo = vc.get_particles_closure(f, geo, mu, rho, p.velocity(), omega=p.omega(), contact_forces=p.contact_forces()) vc.set_coupling_geometry(f, geo) vc.set_coupling_closure(f, clo) # Write output files if i % outf == 0: p.write_mig(outputdir, t) f.write_mig(outputdir, t) # Iterate F = np.zeros(2) # to store inner disk forces mass = np.pi * p.r() ** 2 * rhop time_integration.iterate( f, p, dt, 10, contact_tol=1e-3 * r, external_particles_forces=g * mass, after_sub_iter=lambda n_divide: accumulate(F, n_divide), use_predictor_corrector=False, ) # The FLUID force on the obstacle. NOT get_forces_on_bodies()[0, :]: that # array is per PARTICLE and boundary segments come first, so row 0 is a # boundary segment -- and the fluid populates only the r > 0 grain rows, so # all 242 r == 0 rows are exactly zero by construction. Measured 2026-08-25: # 9100 of 9742 rows nonzero, every one of them a grain, the boundary rows # summing to [0, 0]. boundary_force_by_contributions integrates the traction # over the named boundary and splits it, [p_x, p_y, v_x, v_y]. fb = np.asarray(f.boundary_force_by_contributions("Inner"), float).ravel() Fp, Fv = fb[:2], fb[2:] with open(csv_file, "a") as file1: # columns: total_x; total_y; t; contact_x; contact_y; press_x; press_y; # visc_x; visc_y -- the contact term is kept SEPARATE rather than summed # into the total, so a zero there is visible instead of hidden. Ft = Fp + Fv file1.write(";".join(str(v) for v in (Ft[0], Ft[1], t, F[0], F[1], Fp[0], Fp[1], Fv[0], Fv[1])) + "\n") t += dt i += 1 # %% # Plot # ---- # .. code-block:: shell # # python3 -m migflow.plot.migplot output_2d_drag --actors fluid particles --fluid-field pressure # %% # Artifacts # --------- # - output_2d_drag/animation.mp4