hopper.py 4.74 KB
Newer Older
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
1
# MigFlow - Copyright (C) <2010-2020>
Nathan Coppin's avatar
Nathan Coppin committed
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
# <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 = 1.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 = 100                                    #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
genInitialPosition(outputdir, r, rhop)
Nathan Coppin's avatar
Nathan Coppin committed
73
74
75
p = scontact.ParticleProblem(3,friction_enabled = True,rotation_enabled=True)
p.set_friction_coefficient(0.3,"Sand","Sand")
p.set_friction_coefficient(0.3,"Sand","Steel")
76
p2 = scontact.ParticleProblem(3,friction_enabled = True,rotation_enabled=True)
Nathan Coppin's avatar
Nathan Coppin committed
77
p.read_vtk(outputdir, 0)
78
p2.read_vtk(outputdir, 0)
79
x2 = p2.position()
Nathan Coppin's avatar
Nathan Coppin committed
80
p2.remove_particles_flag((x2[:,1] + p2.r()[:,0] <1.25)*(x2[:,1] + p2.r()[:,0] >1.1)*(x2[:,0] + p2.r()[:,0] <0.18)*(x2[:,0] + p2.r()[:,0] >-0.18))
81
p2.position()[:,1] -= 0.5
82
p.clear_boundaries()
83
p.load_msh_boundaries("mesh" + ".msh", ["Inner","Right","Left","Front","Rear","LockDown"], material = "Steel")
Nathan Coppin's avatar
Nathan Coppin committed
84

Nathan Coppin's avatar
Nathan Coppin committed
85
p.write_vtk(outputdir, 0, 0)
Nathan Coppin's avatar
Nathan Coppin committed
86
87
88
# 
# COMPUTATION LOOP
#
Nathan Coppin's avatar
Nathan Coppin committed
89
opened = 0
Nathan Coppin's avatar
Nathan Coppin committed
90
def accumulate(F):
91
  F+=np.sum(p.get_boundary_forces("Inner"),axis=0)
Nathan Coppin's avatar
Nathan Coppin committed
92
while t < tEnd : 
93
    if t>=10 and opened == 0:
Nathan Coppin's avatar
Nathan Coppin committed
94
95
        opened = 1
        p.clear_boundaries()
96
        p.load_msh_boundaries("mesh" + ".msh", ["Inner","Right","Left","Front","Rear"],material="Steel")
Nathan Coppin's avatar
Nathan Coppin committed
97
#Adding new particles
98
99
    if (max(p.position()[:,1])+max(p.r()) <0.6) and t > 10 and t < 25:
      p.add_particles(p2.position(),p2.r(),p2.mass(),v=p2.velocity(),tag="Sand",forced=p2.forced_flag(),contact_forces=p2.contact_forces())
Nathan Coppin's avatar
Nathan Coppin committed
100
#Solving the contacts
Nathan Coppin's avatar
Nathan Coppin committed
101
    F = np.zeros(3)
102
    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
103
#Removing particles outside the hopper
104
    p.remove_particles_flag( (p.position()[:,1] >-0.6))
Nathan Coppin's avatar
Nathan Coppin committed
105
#Computing mean velocity and writing to file
106
    with open(filename,"a") as file1:
107
        file1.write(str(F[0])+";"+str(F[1])+";"+str(F[2])+";"
Nathan Coppin's avatar
Nathan Coppin committed
108
                    +str(t)+"\n")
Nathan Coppin's avatar
Nathan Coppin committed
109
110
111
    t += dt
    if ii %outf == 0 :
        ioutput = int(ii/outf) + 1
Nathan Coppin's avatar
Nathan Coppin committed
112
        p.write_vtk(outputdir, ioutput, t)      
Nathan Coppin's avatar
Nathan Coppin committed
113
    ii += 1
Nathan Coppin's avatar
Nathan Coppin committed
114
    print("%i : %.2g/%.2g" % (ii, t, tEnd))