hopper.py 5.05 KB
Newer Older
Nathan Coppin's avatar
Nathan Coppin committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
# 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
Nathan Coppin's avatar
Nathan Coppin committed
23
# Study of the drag on a fixed cylinder in a granular flow
Nathan Coppin's avatar
Nathan Coppin committed
24
from migflow import scontact
Nathan Coppin's avatar
Nathan Coppin committed
25
from migflow import fluid
26
from migflow import time_integration
Nathan Coppin's avatar
Nathan Coppin committed
27
28
from migflow import lmgc90Interface
from pylmgc90 import pre
Nathan Coppin's avatar
Nathan Coppin committed
29
30
import numpy as np
import os
Nathan Coppin's avatar
Nathan Coppin committed
31
def genInitialPosition(filename, r, rhop) :
Nathan Coppin's avatar
Nathan Coppin committed
32
33
34
35
36
37
38
39
40
    """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 -- mean radius of the particles
    ly -- domain height
    lx -- domain width 
    rhop -- particles density        
    """
Nathan Coppin's avatar
Nathan Coppin committed
41
    #Particles structure builder
Nathan Coppin's avatar
Nathan Coppin committed
42
    p = scontact.ParticleProblem(3, friction_enabled = True,rotation_enabled=True)
Nathan Coppin's avatar
Nathan Coppin committed
43
    #Load mesh.msh file specifying physical boundaries names
Nathan Coppin's avatar
Nathan Coppin committed
44
    p.load_msh_boundaries("mesh04" + ".msh", ["Inner","Right","Left","Front","Rear","Lockdown"], material = "Steel")
Nathan Coppin's avatar
Nathan Coppin committed
45
    #Particles centre are placed on a regular grid
Nathan Coppin's avatar
Nathan Coppin committed
46
    # Add a grain at each centre position
Nathan Coppin's avatar
Nathan Coppin committed
47
48
    for y in np.arange(1+r, 2.7, 2*r*1.05):
            for x in np.arange(-0.19+r*1.05, 0.19-r*1.05, 2*r*1.05):
Nathan Coppin's avatar
Nathan Coppin committed
49
50
51
                for z in np.arange(0+r*1.05,0.05-1.05*r,2*r*1.05):
                    r1 = np.random.uniform(.95*r, 1.05*r)
                    p.add_particle((x, y,z), r1, (4/3)*r1**3 * np.pi * rhop,"Sand");  
Nathan Coppin's avatar
Nathan Coppin committed
52
    p.write_vtk(outputdir, 0, 0)
Nathan Coppin's avatar
Nathan Coppin committed
53
54
# Physical parameters
g = np.array([0,-9.81,0])
Nathan Coppin's avatar
Nathan Coppin committed
55
r = 1e-2
Nathan Coppin's avatar
Nathan Coppin committed
56
57
58
59
rhop = 2500
# Numerical parameters
dt = 1e-3                                      #time step
t = 0                                          #initial time
60
tEnd = 40                                      #final time
Nathan Coppin's avatar
Nathan Coppin committed
61
ii = 0                                         #iteration number
Nathan Coppin's avatar
Nathan Coppin committed
62
outf = 1000                                    #iterations between data frames  
Nathan Coppin's avatar
Nathan Coppin committed
63
# Define output directory
Nathan Coppin's avatar
Nathan Coppin committed
64
outputdir = "output"
Nathan Coppin's avatar
Nathan Coppin committed
65
66
if not os.path.isdir(outputdir) :
    os.makedirs(outputdir)
Nathan Coppin's avatar
Nathan Coppin committed
67
#Creating file for drag data
Nathan Coppin's avatar
Nathan Coppin committed
68
filename = outputdir + "/Drag.csv"
Nathan Coppin's avatar
Nathan Coppin committed
69
70
71
# 
# PARTICLE PROBLEM
#
Nathan Coppin's avatar
Nathan Coppin committed
72
fric = friction[num][0]
Nathan Coppin's avatar
Nathan Coppin committed
73
genInitialPosition(outputdir, r, rhop)
Nathan Coppin's avatar
Nathan Coppin committed
74
75
76
77
78
79
if use_lmgc :
  friction = 0.3                                #friction coefficient
  lmgc90Interface.scontactTolmgc90(outputdir, 3, 0, fric)
  p = lmgc90Interface.ParticleProblem(3)
else : 
  p = scontact.ParticleProblem(3,friction_enabled = True,rotation_enabled=True)
Nathan Coppin's avatar
Nathan Coppin committed
80
81
  p.set_friction_coefficient(0.3,"Sand","Sand")
  p.set_friction_coefficient(0.3,"Sand","Steel")
82
p2 = scontact.ParticleProblem(3,friction_enabled = True,rotation_enabled=True)
Nathan Coppin's avatar
Nathan Coppin committed
83
p.read_vtk(outputdir, 0)
84
p2.read_vtk(outputdir, 0)
Nathan Coppin's avatar
Nathan Coppin committed
85
86
87
88
89
90
91
state = p.state()
state2 = p2.state()
p2.remove_particles_flag((state2.x[:,1] + p2.r()[:,0] <1.25))
p2.remove_particles_flag((state2.x[:,1] + p2.r()[:,0] >1.1))
p2.remove_particles_flag((state2.x[:,0] + p2.r()[:,0] <0.18))
p2.remove_particles_flag((state2.x[:,0] + p2.r()[:,0] >-0.18))
state2.x[:,1] -= 0.5
92
p.clear_boundaries()
93
p.load_msh_boundaries("mesh" + ".msh", ["Inner","Right","Left","Front","Rear","LockDown"], material = "Steel")
Nathan Coppin's avatar
Nathan Coppin committed
94

Nathan Coppin's avatar
Nathan Coppin committed
95
p.write_vtk(outputdir, 0, 0)
Nathan Coppin's avatar
Nathan Coppin committed
96
97
98
# 
# COMPUTATION LOOP
#
Nathan Coppin's avatar
Nathan Coppin committed
99
opened = 0
Nathan Coppin's avatar
Nathan Coppin committed
100
def accumulate(F):
Nathan Coppin's avatar
Nathan Coppin committed
101
  F+=np.sum(p.get_boundary_impulses("Inner"),axis=0)
Nathan Coppin's avatar
Nathan Coppin committed
102
while t < tEnd : 
103
    if t>=10 and opened == 0:
Nathan Coppin's avatar
Nathan Coppin committed
104
105
        opened = 1
        p.clear_boundaries()
106
        p.load_msh_boundaries("mesh" + ".msh", ["Inner","Right","Left","Front","Rear"],material="Steel")
Nathan Coppin's avatar
Nathan Coppin committed
107
#Adding new particles
Nathan Coppin's avatar
Nathan Coppin committed
108
109
    if (max(state.x[:,1])+max(p.r()) <0.6) and t > 10 and t < 25:
      p.add_particles(state2.x,p2.r(),p2.mass(),v=state2.v,tag="Sand",forced=p2.forced_flag(),contact_forces=p2.contact_forces())
Nathan Coppin's avatar
Nathan Coppin committed
110
#Solving the contacts
Nathan Coppin's avatar
Nathan Coppin committed
111
    F = np.zeros(3)
112
    time_integration.iterate(None,p,dt,min_nsub=1,contact_tol=1e-6,external_particles_forces=g*p.mass(),after_sub_iter= lambda subdt : accumulate(F))
Nathan Coppin's avatar
Nathan Coppin committed
113
#Removing particles outside the hopper
Nathan Coppin's avatar
Nathan Coppin committed
114
    p.remove_particles_flag( (state.x[:,1] >-0.6))
Nathan Coppin's avatar
Nathan Coppin committed
115
#Computing mean velocity and writing to file
116
    with open(filename,"a") as file1:
117
        file1.write(str(F[0])+";"+str(F[1])+";"+str(F[2])+";"
Nathan Coppin's avatar
Nathan Coppin committed
118
                    +str(t)+"\n")
Nathan Coppin's avatar
Nathan Coppin committed
119
120
121
    t += dt
    if ii %outf == 0 :
        ioutput = int(ii/outf) + 1
Nathan Coppin's avatar
Nathan Coppin committed
122
        p.write_vtk(outputdir, ioutput, t)      
Nathan Coppin's avatar
Nathan Coppin committed
123
    ii += 1
Nathan Coppin's avatar
Nathan Coppin committed
124
    print("%i : %.2g/%.2g" % (ii, t, tEnd))