melangeur.py

# MigFlow - Copyright (C) <2010-2018>
# <Universite catholique de Louvain (UCL), Belgium
#  Universite de Montpellier, France>
# 	
# List of the contributors to the development of MigFlow: see AUTHORS file.
# Description and complete License: see LICENSE file.
# 	
# This program (MigFlow) is free software: 
# you can redistribute it and/or modify it under the terms of the GNU Lesser General 
# Public License as published by the Free Software Foundation, either version
# 3 of the License, or (at your option) any later version.
# 
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU Lesser General Public License for more details.
# 
# You should have received a copy of the GNU Lesser General Public License
# along with this program (see COPYING and COPYING.LESSER files).  If not, 
# see <http://www.gnu.org/licenses/>.

#!/usr/bin/env python

# TESTCASE DESCRIPTION
# Grains in circular shear flow.
# This example shows how to set boundary conditions as a function of some parametres.
# A boolean parametre give the possibility to use lmgc90 to solve contacts or not.

from migflow import fluid
from migflow import scontact

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

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

def genInitialPosition(filename, r, rout, rin, rhop) :
    """Create a container to initialize particles and write particles in an initial output file.

    Keyword arguments:
    filename -- directory to write the output file
    r -- max particle radius
    rout -- outer radius of the drum
    rhop -- particle density
    """
    # Grains structure builder
    p = scontact.ParticleProblem(2)
    # Load mesh.msh file specifying physical boundaries names
    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]
    # Add a grain at each centre position
    for i in range(x.shape[0]) :
        if y[i]<0:
            rhop1 = random.choice([rhop*.9,1.1*rhop,rhop])
            p.add_particle((x[i], y[i]), r, r**2 * np.pi * rhop1)
    p.write_vtk(filename,0,0)
 
# Physical parameters
g = 0                                          # gravity
rho = 1.253e3                                  # fluid density
rhop = 1000                                    # grains density
nu = 1e-4                                      # kinematic viscosity

# Numerical parameters
dt = 1e-2                                      # time step
tEnd = 50                                      # final time

# Geometry parameters
outf = 5                                       # number of iterations between output files
rout = 0.0254                                  # outer radius
rin = 0.0064                                   # inner radius
r = 397e-6/2                                   # grains radius

#
# PARTICLE PROBLEM
#
# Object particles creation
genInitialPosition(outputdir, r, rout, rin, rhop)
p = scontact.ParticleProblem(2)
p.read_vtk(outputdir,0)

# Initial time and iteration
t = 0
ii = 0

#
# FLUID PROBLEM
#
fluid = fluid.FluidProblem(2,g,[nu*rho],[rho])
# Set the mesh geometry for the fluid computation
fluid.load_msh("mesh.msh")
fluid.set_strong_boundary("Outer",0,lambda x : -x[:, 1])
fluid.set_strong_boundary("Outer",1,lambda x : x[:, 0])
fluid.set_strong_boundary("Inner",0,0)
fluid.set_strong_boundary("Inner",1,0)
fluid.set_strong_boundary("PtFix",2,0)
fluid.set_weak_boundary("Inner","Wall")
fluid.set_weak_boundary("Outer","Wall")
# Set location of the grains in the mesh and compute the porosity in each computation cell
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.export_vtk(outputdir,0,0)

tic = time.time()
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())

# 
# COMPUTATION LOOP
#
while t < tEnd :
    # Fluid solver
    fluid.implicit_euler(dt)

    # Computation of the new velocities
    forces = fluid.compute_node_force(dt)
    vn = p.velocity() + forces*dt/p.mass()
    vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
    # Number of contact 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()))
    # NLGS iterations
    for i in range(nsub) :
        tol = 1e-6                                # Numerical tolerance for grains intersection
        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.time()-tic))