melangeur.py

#!/usr/bin/env python
from marblesbag import fluid as fluid
from marblesbag import scontact2

import numpy as np
import os
import time
import shutil
import random

outputdir = "output"
if not os.path.isdir(outputdir) :
    os.makedirs(outputdir)

#grains creation and placement + physical boundaries definition
'''r is the radius of the grains
rout is the outer radius of the geometry
rin is the inner radius of the geometry
rhop is the particles density'''
def genInitialPosition(filename, r, rout, rin, rhop) :
    p = scontact2.ParticleProblem()
    #Loading of the mesh.msh file specifying physical boundaries name
    p.load_msh_boundaries("mesh.msh", ["Outer", "Inner"])
    #Definition of the points where the grains are located
    x = np.arange(rout, -rout, 2 * -r)
    x, y = np.meshgrid(x, x)
    R2 = x**2 + y**2
    #condition to be inside the outer boundary
    keep = R2 < (rout - r)**2
    x = x[keep]
    y = y[keep]
    R2 = x**2 + y**2
    #condition to be outside the inner boundary
    keep = R2 > (rin + r)**2
    x = x[keep]
    y = y[keep]
    for i in range(x.shape[0]) :
        if y[i]<0:
            rhop1 = random.choice([rhop*.9,1.1*rhop,rhop])
            #Addition of an particle object at each point
            p.add_particle((x[i], y[i]), r, r**2 * np.pi * rhop1);
    p.write_vtk(filename,0,0)


t = 0
ii = 0

#physical parameters
g =  0                                  # gravity
rho = 1.253e3                           # fluid density
rhop = 1000                             # grains density
nu = 1e-4                               # kinematic viscosity
tEnd = 50                               # final time

#numerical parameters
h = 0.002                               # mesh size
dt = 1e-2                               # time step
epsilon = 1e-5                          # stabilisation parametre

#geometry parameters
rout = 0.0254                           # outer radius
rin = 0.0064                            # inner radius
r = 397e-6/2                            # grains radius
shutil.copy("mesh.msh", outputdir +"/mesh.msh")

#Object particles creation
genInitialPosition(outputdir, r, rout, rin, rhop)
p = scontact2.ParticleProblem()
p.read_vtk(outputdir,0)

outf = 5                                # number of iterations between output files

#Object fluid creation + Boundary condition of the fluid (field 0 is horizontal velocity; field 1 is vertical velocity; field 2 is pressure)
#Format: strong_boundaries = [(Boundary tag, Fluid field, Value)
strong_boundaries = [("Outer",0,lambda x : -x[:, 1]),("Outer",1,lambda x : x[:, 0]),("Inner",0,0.),("Inner",1,0.),("PtFix",2,0.)]
fluid = fluid.fluid_problem("mesh.msh",g,nu*rho,rho,epsilon,strong_boundaries)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.export(outputdir,0,0)

tic = time.clock()
#Computation loop
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
while t < tEnd :
    #Fluid solver
    fluid.implicit_euler(dt)
    forces = fluid.compute_node_force(dt)
    #Computation of the new velocities
    vn = p.velocity() + forces * dt / p.mass()
    vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
    #number of sub time step
    nsub = max(1, int(np.ceil((vmax * dt * 4)/min(p.r()))))
    print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
    #Contact solver
    for i in range(nsub) :
        tol = 1e-6
        p.iterate(dt/nsub, forces, tol)
    #Localisation of the grains in the fluid
    fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
    t += dt
    #Output files writting
    if ii %outf == 0 :
        ioutput = int(ii/outf) + 1
        p.write_vtk(outputdir, ioutput, t)
        fluid.export_vtk(outputdir, t, ioutput)
    ii += 1
    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.clock() - tic))