hopper.py 5.09 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
25
26
27
28
from migflow import scontact
from migflow import fluid
from migflow import time_integration
import numpy as np
import os
Nathan Coppin's avatar
Nathan Coppin committed
29
30
31
32
33
34
35
36
37
38
39
40
41
def genInitialPosition(filename, r, 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 -- mean radius of the particles
    ly -- domain height
    lx -- domain width 
    rhop -- particles density        
    """
    #Particles structure builder
    p = scontact.ParticleProblem(2, friction_enabled = True,rotation_enabled=True)
    #Load mesh.msh file specifying physical boundaries names
Nathan Coppin's avatar
Nathan Coppin committed
42
43
44
    p.load_msh_boundaries("mesh.msh", ["Inner"], material = "Steel")
    p.load_msh_boundaries("mesh.msh", ["Right","Left","RightUp","RightDown","LeftUp","LeftDown"], material = "Plexi")
    p.load_msh_boundaries("mesh.msh", ["LockDown"], material = "Steel")
Nathan Coppin's avatar
Nathan Coppin committed
45
46
47
48
49
50
51
    #Particles centre are placed on a regular grid
    # Add a grain at each centre position
    for y in np.arange(1+r, 2, 2*r*1.05):
            for x in np.arange(-0.295+r*1.05, 0.295-r*1.05, 2*r*1.05):
                r1 = np.random.uniform(.95*r, 1.05*r)
                p.add_particle((x, y), r1, r1**2 * np.pi * rhop,"Sand");  
    p.write_vtk(outputdir, 0, 0)
Nathan Coppin's avatar
Nathan Coppin committed
52
53
# Physical parameters
g = np.array([0,-9.81])
Nathan Coppin's avatar
Nathan Coppin committed
54
r = 8e-3
Nathan Coppin's avatar
Nathan Coppin committed
55
rhop = 2500
56
mupart= 0.12
Nathan Coppin's avatar
Nathan Coppin committed
57
58
59
# Numerical parameters
dt = 1e-3                                      #time step
t = 0                                          #initial time
60
tEnd = 13                                      #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
67
68
69
70
71
if not os.path.isdir(outputdir) :
    os.makedirs(outputdir)
#Creating file for drag data
filename = outputdir + "/Drag.csv"
# 
# PARTICLE PROBLEM
#
Nathan Coppin's avatar
Nathan Coppin committed
72
73
74
75
76
genInitialPosition(outputdir, r, rhop)
p = scontact.ParticleProblem(2,friction_enabled = True,rotation_enabled=True)
p2 = scontact.ParticleProblem(2,friction_enabled = True,rotation_enabled=True)
p.read_vtk(outputdir, 0)
p2.read_vtk(outputdir, 0)
Nathan Coppin's avatar
Nathan Coppin committed
77
78
79
80
81
82
83
state = p.state()
state2 = p2.state()
p2.remove_particles_flag((state2.x[:,1] + p2.r()[:,0] <1.15))
p2.remove_particles_flag((state2.x[:,1] + p2.r()[:,0] >1.1))
p2.remove_particles_flag((state2.x[:,0] + p2.r()[:,0] <0.15))
p2.remove_particles_flag((state2.x[:,0] + p2.r()[:,0] >-0.15))
state2.x[:,1] -= 0.6
Nathan Coppin's avatar
Nathan Coppin committed
84
p.clear_boundaries()
Nathan Coppin's avatar
Nathan Coppin committed
85
86
87
p.load_msh_boundaries("mesh.msh", ["Inner"], material = "Steel")
p.load_msh_boundaries("mesh.msh", ["Right","Left","RightUp","RightDown","LeftUp","LeftDown"], material = "Plexi")
p.load_msh_boundaries("mesh.msh", ["LockDown"], material = "Steel")
88
89
90
p.set_friction_coefficient(mupart,"Sand","Sand")
p.set_friction_coefficient(0.1,"Sand","Steel")
p.set_friction_coefficient(0.1,"Sand","Plexi")
Nathan Coppin's avatar
Nathan Coppin committed
91
92
93
94
p.write_vtk(outputdir, 0, 0)
# 
# COMPUTATION LOOP
#
Nathan Coppin's avatar
Nathan Coppin committed
95
96
opened = 0
def accumulate(F):
Nathan Coppin's avatar
Nathan Coppin committed
97
  F+=np.sum(p.get_boundary_forces("Inner"),axis=0)
Nathan Coppin's avatar
Nathan Coppin committed
98
while t < tEnd : 
99
    if t>=3 and opened == 0:
Nathan Coppin's avatar
Nathan Coppin committed
100
101
        opened = 1
        p.clear_boundaries()
Nathan Coppin's avatar
Nathan Coppin committed
102
103
        p.load_msh_boundaries("mesh.msh", ["Inner"],material="Steel")
        p.load_msh_boundaries("mesh.msh", ["Right","Left","RightUp","RightDown","LeftUp","LeftDown"],material="Plexi")
Nathan Coppin's avatar
Nathan Coppin committed
104
#Adding new particles
Nathan Coppin's avatar
Nathan Coppin committed
105
106
    if (max(state.x[:,1])+max(p.r()) <0.5) and t > 3 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
107
#Solving the contacts
Nathan Coppin's avatar
Nathan Coppin committed
108
    F = np.zeros(2)
109
    time_integration.iterate(None,p,dt,min_nsub=1,contact_tol=1e-7,external_particles_forces=g*p.mass(),after_sub_iter= lambda subdt : accumulate(F))
Nathan Coppin's avatar
Nathan Coppin committed
110
#Removing particles outside the hopper
Nathan Coppin's avatar
Nathan Coppin committed
111
    p.remove_particles_flag( (state.x[:,1] >-0.6))
Nathan Coppin's avatar
Nathan Coppin committed
112
113
#Computing mean velocity and writing to file
    with open(filename,"a") as file1:
114
        file1.write(str(F[0])+";"+str(F[1])+";"
Nathan Coppin's avatar
Nathan Coppin committed
115
                    +str(t)+"\n")
Nathan Coppin's avatar
Nathan Coppin committed
116
117
118
119
120
121
    t += dt
    if ii %outf == 0 :
        ioutput = int(ii/outf) + 1
        p.write_vtk(outputdir, ioutput, t)      
    ii += 1
    print("%i : %.2g/%.2g" % (ii, t, tEnd))