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.1                                     #stream velocity
rhop = 2500                                    #particles density
rho = 1000                                     #fluid density
nu = 1e-6                                      #fluid kinematic viscosity
g = np.array([0,-9.81])                        #gravity
# Numerical parameters
dt = 1e-3                                      #time step
t = 0                                          #initial time
tEnd = 13                                      #final time
ii = 0                                         #iteration number
outf = 10                                      #iterations between data frames  
# Define output directory
outputdir = "output" 
if not os.path.isdir(outputdir) :
    os.makedirs(outputdir)
filename = outputdir + "/Drag.csv"
# 
# PARTICLE PROBLEM
#
#Injected particles
p = scontact.ParticleProblem(2,True,True)
p.read_vtk("depot/output",10)
state = p.state()
p.remove_particles_flag((state.x[:,1]-p.r()[:,0]>0.5)*(state.x[:,1]+p.r()[:,0] <0.53))
#Flow particles
p2 = scontact.ParticleProblem(2,True,True)
p2.read_vtk("depot/output",10)
state2 = p2.state()
p2.remove_particles_flag(state2.x()[:,1]+p2.r()[:,0] <0.5)
p2.clear_boundaries()
p2.load_msh_boundaries("mesh.msh", ["Inner","Left","Right"],material = "Steel") 
p2.set_friction_coefficient(0.3,"Sand","Sand")
p2.set_friction_coefficient(0.3,"Sand","Steel")
p2.write_vtk(outputdir, 0, 0)
#
# FLUID PROBLEM
#
fluid = fluid.FluidProblem(2,g,[nu*rho],[rho],drag_in_stab=0)
fluid.load_msh("mesh.msh")
fluid.set_wall_boundary("Inner")
fluid.set_wall_boundary("Left", velocity = [0,0])
fluid.set_wall_boundary("Right", velocity = [0,0])
fluid.set_open_boundary("Bottom", velocity = [0,vit])
fluid.set_open_boundary("Top",pressure = 0)
fluid.set_particles(p2.mass(), p2.volume(), state2, p2.contact_forces(),init=True)
fluid.export_vtk(outputdir, 0.,0)
# 
# COMPUTATION LOOP
#
def accumulate(F):
  F+=np.sum(p2.get_boundary_impulses("Inner"),axis=0)
while t < tEnd : 
     p2.forced_flag()[(state2.x[:,1] + p2.r()[:,0] <= -0.5)] = 1
     state2.v[p2.forced_flag()==1,:] = [0, vit]
     state2.omega[p2.forced_flag()==1] = 0
     if  max(state2.x[:,1]+p2.r()[:,0]) < 0.5 and t < 25 :
       p2.add_particles(state.x,p.r(),p.mass(),v=state.v,tag="Sand",forced=p.forced_flag(),contact_forces=p.contact_forces())
     p2.remove_particles_flag( (state2.x[:,1] + p2.r()[:,0] >-0.55))
     F = np.zeros(2)
     time_integration.particle_changed(fluid,p2)
     time_integration.iterate(fluid,p2,dt,1,contact_tol=5e-7,external_particles_forces=g*p2.mass(),after_sub_iter=lambda subdt : accumulate(F))
     with open(filename,"a") as file1:
        F += fluid.boundary_force()[0,:]
        file1.write( str(F[0])+";"+str(F[1])+";"+str(t)+"\n")
     t += dt   
     # Output files writing
     if ii %outf == 0 :
         ioutput = int(ii/outf) + 1
         p2.write_vtk(outputdir, ioutput, t)
         fluid.export_vtk(outputdir+"/Flow", t, ioutput)       
     ii += 1
     print("%i : %.2g/%.2g" % (ii, t, tEnd))