drop.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
# This tescase presents the fall of a cloud made of particles in a viscous fluid (Stokes cloud).
# Physical parameters for the drops are the ones presented by Metzger et al. (2007) 
# "Falling clouds of particles in viscous fluids"

from migflow import fluid
from migflow import scontact
from migflow import time_integration

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

def genInitialPosition(filename, r, rout, rhop, compacity) :
    """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 -- max radius of the particles
    rout -- outer radius of the cloud
    rhop -- particles density        
    compacity -- initial compacity in the cloud
    """
    # Particles structure builder
    p = scontact.ParticleProblem(2)
    # Load mesh.msh file specifying physical boundaries names
    p.load_msh_boundaries("mesh.msh", ["Top", "Bottom","Lateral"])

    # Space between the particles to obtain the expected compacity
    N = compacity*4*rout**2/(np.pi*r**2)
    e = 2*rout/(N)**(1./2.)

    # Definition of the points where the particles are located
    x = np.arange(rout, -rout, -e)
    x, y = np.meshgrid(x, x)
    R2 = x**2 + y**2
    keep = R2 < (rout - r)**2
    x = x[keep]
    y = y[keep]

    for i in range(x.shape[0]) :
        # Add particles object with properties defined above
        p.add_particle((x[i]+2*random.uniform(-e/2+r,e/2-r), y[i]+2*random.uniform(-e/2+r,e/2-r)), r, r**2 * np.pi * rhop)
    
    # Vertical shift of the particles to the top of the box
    p.position()[:, 1] += 0.52
    p.write_vtk(filename,0,0)

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

# Physical parameters
g = -9.81                                       # gravity
rhop = 2450                                     # particles density
r = 154e-6/2.5                                  # particles radii
compacity = 0.20                                # solid volume fraction in the drop
rho = 1030                                      # fluid density
nu = 1.17/rho                                   # kinematic viscosity
rout = 3.3e-3                                   # cloud radius
mu = nu*rho                                     # dynamic viscosity
print("RHOP = %g, r = %g, RHO = %g" % (rhop,r,rho))

# Numerical parameters
outf = 1                                        # number of iterations between output files
dt = 2e-2                                       # time step
tEnd = 100                                      # final time

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

# Initial time and iteration
t = 0
ii = 0
jj = 0
tEnd1 = 0.2

#
# FLUID PROBLEM
#
fluid = fluid.FluidProblem(2,g,[nu*rho],[rho],drag_in_stab=1)
# Set the mesh geometry for the fluid computation
fluid.load_msh("mesh.msh")
fluid.set_wall_boundary("Top", pressure=0)
fluid.set_wall_boundary("Bottom")
fluid.set_wall_boundary("Lateral")
# Set location of the particles 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(),init=True)

t = 0
fluid.solution()[:,2] = (fluid.coordinates()[:,1]-0.6)*rho*g
fluid.export_vtk(outputdir,0,0)
tic = time.time()

G = np.zeros_like(p.velocity())
G[:,1] = p.mass()[:,0]*g
#
# COMPUTATION LOOP
#
while t < tEnd :

    if ((ii+1)%15==0):
        number_p = fluid.n_particles
        position_p = fluid.particle_position()
        volume_p = fluid.particle_volume()
    if (ii%15==0 and ii != 0):
        fluid.adapt_mesh(2e-2,1e-3,10000,old_n_particles=number_p,old_particle_position=position_p,old_particle_volume=volume_p)

    time_integration.iterate(fluid, p, dt, 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
    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))