drag.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
# Study of the drag on a fixed disk in a 2D immersed granular flow
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

# Physical parameters
vit = -0.48                                    #stream velocity
rhop = 2500                                    #particles density
rho = 500                                     #fluid density
nu = 1e-8                                      #fluid kinematic viscosity
g = np.array([0,-21.81])
# Numerical parameters
dt = 1e-3                                      #time step
t = 0                                          #initial time
tEnd = 13                                      #final time
ii = 0                                         #iteration number
outf = 1                                      #iterations between data frames  

# Define output directory
outputdir = "ellipseForced"
if not os.path.isdir(outputdir) :
    os.makedirs(outputdir)
filename = outputdir + "/Drag.csv"
if os.path.exists(filename):
  print("Remove old results? (y/n)\n")
  decision = input()
  if decision.lower() == 'y':
    os.remove(filename)

# 
# PARTICLE PROBLEM
#

#Injected particles
p = scontact.ParticleProblem(2,True,True)
p.read_vtk("depot/depotEllipse",450) 
p.remove_particles_flag( (p.position()[:,1]  < -0.2339))                   
p.position()[:,:] = p.position()[:,(1,0)]
p.position()[:,:] += [0.38005,0.647]
p.remove_particles_flag( (p.position()[:,1]  > 0.507))
p.velocity()[:,:] = [0,0] 
flipped = np.copy(p.position())
flipped[:,0] = -flipped[:,0]
p.contact_forces()[:,:] = 0
p.omega()[:] = 0
p.forced_flag()[:] = 1
p.mass()[:] *= 0.05
p.write_vtk(outputdir + "/In_particles", 0, 0)
#Injected particles


#Flow particles
p2 = scontact.ParticleProblem(2,True,True)
p2.read_vtk(outputdir + "/In_particles",0)
p2.clear_boundaries()
p2.load_msh_boundaries("mesh.msh", ["Inner"],material = "Steel") 
p2.load_msh_boundaries("mesh.msh", ["Left", "Right"],material = "Plexi")
#p2.remove_particles_flag( (p2.position()[:,1]  < 0.5-2*p2.r()[:,0]))
#p2.add_particles(p.position(),p.r(),p.mass(),v=p.velocity(),tag="Sand",forced=p.forced_flag(),contact_forces=p.contact_forces())
p2.add_particles(p.position()[:,:] - [0,0.291],p.r(),p.mass(),v=np.zeros_like(p.velocity()),tag="Sand",forced=np.zeros_like(p.forced_flag()),contact_forces=p.contact_forces())
p2.add_particles(flipped - [0,0.291*2],p.r(),p.mass(),v=np.zeros_like(p.velocity()),tag="Sand",forced=np.zeros_like(p.forced_flag()),contact_forces=p.contact_forces())
p2.add_particles(p.position()[:,:] - [0,0.291*3],p.r(),p.mass(),v=np.zeros_like(p.velocity()),tag="Sand",forced=np.zeros_like(p.forced_flag()),contact_forces=p.contact_forces())
p2.add_particles(flipped - [0,0.291*4],p.r(),p.mass(),v=np.zeros_like(p.velocity()),tag="Sand",forced=np.zeros_like(p.forced_flag()),contact_forces=p.contact_forces())
p2.remove_particles_flag( (p2.position()[:,1]**2 + p2.position()[:,0]**2  > (0.01905 + max(p2.r()))**2))
p2.r()[:] *= 0.99


p2.set_friction_coefficient(0.28,"Sand","Sand")
p2.set_friction_coefficient(0.28,"Sand","Steel")
p2.set_friction_coefficient(0.14,"Sand","Plexi")
p2.set_tangent_boundary_velocity("Left",0)#vit)
p2.set_tangent_boundary_velocity("Right",0)#-vit)
p2.forced_flag()[p2.position()[:,1] + 2*p2.r()[:,0] >= 0.5] = 1
p2.contact_forces()[:,:] = 0
p2.velocity()[:,0] = 0
p2.velocity()[p2.position()[:,1] < 0,1] = 0
p2.omega()[:] = 0
p2.mass()[:,0] *= 0.05
p2.set_get_momentum(1) 
p2.write_vtk(outputdir+"/Flow", 0, 0)

#
# FLUID PROBLEM
#

#fluid = fluid.FluidProblem(2,0,[nu*rho],[rho])
#Set the mesh geometry for the fluid computation
#fluid.load_msh("mesh.msh")
#fluid.set_wall_boundary("Inner")
#fluid.set_wall_boundary("Left", velocity = [0,vit])
#fluid.set_wall_boundary("Right", velocity = [0,vit])
#fluid.set_open_boundary("Bottom", velocity = [0,vit])
#fluid.set_open_boundary("Top",pressure = 0)
#fluid.set_particles(p2.mass(), p2.volume(), p2.position(), p2.velocity(), p2.contact_forces(),init=True)
#fluid.export_vtk(outputdir+"/Flow", 0.,0)

# 
# COMPUTATION LOOP
#++	
print(p2.n_particles())
tic = time.time()
vitC = -0.01
while t < tEnd : 
     momentum = np.zeros_like(p2.momentum())
     if vitC >= vit and t>=500*dt:
       vitC += 2*dt*vit
     vitB = vitC#*0.72 if int((t+dt)/dt)%2==0 else vitC#-0.002 if int((t/dt))%3==0 else vitC+0.002
     vitC2 = vitC#*72.0 if int(t/dt)%2==0 else vitC
     print("vitC", vitC, "vitB", vitB)
     p2.forced_flag()[p2.position()[:,1] + 2*p2.r()[:,0] < 0.5] = 0
     p2.forced_flag()[(p2.position()[:,1] + 2*p2.r()[:,0] <= -0.5)+(p2.position()[:,1] + 2*p2.r()[:,0] >= 0.5)] = 1
     p2.velocity()[p2.forced_flag()==1,:] = [0, vitC2]
     p2.velocity()[(p2.position()[:,1] + 2*p2.r()[:,0] <= -0.5),1] = vitB
     p2.omega()[p2.forced_flag()==1] = 0
     if  max(p2.position()[:,1]) < 0.5 :
       p2.add_particles(p.position() if random.choice([True,False]) else flipped,p.r(),p.mass(),v=p.velocity(),tag="Sand",forced=p.forced_flag(),contact_forces=p.contact_forces())
     p2.remove_particles_flag( (p2.position()[:,1] + p2.r()[:,0] >-0.551))
     #time_integration.particle_changed(fluid,p2)
     G = np.zeros_like(p2.velocity())
     G[:,0] = p2.mass()[:,0]*g[1]
     time_integration.iterate(None,p2,dt,20,contact_tol=3e-8,momentum=momentum)
     #momentum += fluid.boundary_force()[0:3,:]
     print(momentum)
     speedAbove = np.mean(p.velocity()[(p.position()[:,1] < 0.45)*(p.position()[:,1] > 0.3),1])
     speedBelow = np.mean(p.velocity()[(p.position()[:,1] > -0.45)*(p.position()[:,1] < -0.3),1])
     with open(filename,"a") as file1:
        file1.write(str(momentum[0,0]) + ";" + str(momentum[0,1]) + ";" + str(momentum[1,0]) + ";" 
                    + str(momentum[1,1]) + ";" + str(momentum[2,0]) + ";" + str(momentum[2,1]) + ";" 
                    + str(t) + ";" + str(0 if np.isnan(speedAbove) else speedAbove) + ";" + str(0 if np.isnan(speedBelow) else speedBelow) + "\n")
     t += dt   
     # Output files writing
     if ii %outf == 0 :
         ioutput = int(ii/outf) + 1
         p2.write_vtk(outputdir+"/Flow", ioutput, t)
         #fluid.export_vtk(outputdir+"/Flow", t, ioutput)       
     ii += 1
     print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))