mill.py

# Migflow - Copyright (C) <2010-2018>
# <Universite catholique de Louvain (UCL), Belgium
#  Universite de Montpellier, France>
# 	
# List of the contributors to the development of Marblesbag: 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
from migflow import scontact
from migflow import fluid
from migflow import time_integration

import numpy as np
import os
import time
import shutil
import random
outputdir = "output"
if not os.path.isdir(outputdir) :
    os.makedirs(outputdir)


def genInitialPosition(filename, r, rout, rhop) :
    """Set all the particles centre positions and create the particles objects to add in the computing structure

    Keyword arguments:
    filename -- name of the output file
    r -- radius of the particles
    rout -- the outer radius of the geometry
    rin -- the inner radius of the geometry
    rhop -- the particles density
    """
    p = scontact.ParticleProblem(2)
    #Loading of the mesh.msh file specifying physical boundaries name and material
    p.load_msh_boundaries("mesh.msh", ["Outer","Top",],material="PVC")
    #Definition of the points where the grains are located
    x = np.arange(rout-1.2*r, -rout+1.2*r, -2*r*1.05)
    x, y = np.meshgrid(x, x)
    R2 = x**2 + y**2
    #condition to be inside the outer boundary
    keep = R2 < (rout - r*1.05)**2
    x = x[keep]
    y = y[keep]
    R2 = x**2 + y**2
    for i in range(x.shape[0]) :
        if y[i]<0.005:
            r1 = np.random.normal(r,0.05*r)
            #Addition of an particle object at each point
            p.add_particle((x[i], y[i]), r1, r1**2 * np.pi * rhop,"Glass");
    p.write_vtk(filename,0,0)

t = 0
ii = 0
#physical parameters
g    = -9.81                            # gravity
rhop = 2500                             # grains density
rho = 1000
nu = 1e-6
tEnd = 10                               # final time

#numerical parameters
dt = 1e-3                               # time step

#geometry parameters
rout = 0.05                             # outer radius
r  = 1e-3                               # grains radius
Fr = 8e-3				                # Froude number
v = 0.0715                              # Corresponding tangent velocity

shutil.copy("mesh.msh", outputdir +"/mesh.msh")

#
# PARTICLE PROBLEM
#
# Initialise particles
genInitialPosition(outputdir, r, rout, rhop)
p = scontact.ParticleProblem(2,True)
p.read_vtk(outputdir,0)
outf = 250                                       # number of iterations between output files
#Taking friction into account
p.set_friction_coefficient(.4,"Glass","Glass")   # between particles
p.set_friction_coefficient(.5,"Glass","PVC")     # between particles and boundaries

p.set_tangent_boundary_velocity("Outer",-v)      # setting the tangent velocity of the rotating mill
p.set_tangent_boundary_velocity("Top",-v)        # setting the tangent velocity of the rotating mill

#
# 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]/rout)*-v )
fluid.set_strong_boundary("Outer",1,lambda x : (-x[:,0]/rout)*-v)
fluid.set_wall_boundary("Top",0)
# 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(),p.contact_forces())
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces())
fluid.solution()[:,2] = g*rho*(fluid.coordinates()[:,1] - rout/2)
fluid.export_vtk(outputdir, 0.,0)

G = np.zeros_like(p.velocity())
G[:,1] = p.mass()[:,0]*g

#
# COMPUTATION LOOP
#
while t < tEnd :
    time_integration.iterate(fluid, p, dt, min_nsub=15, external_particles_forces=G)
    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