inject2fGrains.py 4.49 KB
Newer Older
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
1
# MigFlow - Copyright (C) <2010-2020>
Matthieu Constant's avatar
Matthieu Constant committed
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
# <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
from migflow import fluid
from migflow import scontact
from migflow import lmgc90Interface
from pylmgc90 import pre
import numpy as np
import os
import time
import shutil
import random

Matthieu Constant's avatar
Matthieu Constant committed
33
outputdir = "output1"
Matthieu Constant's avatar
Matthieu Constant committed
34
35
36
37
38
if not os.path.isdir(outputdir) :
    os.makedirs(outputdir)

t = 0
ii = 0
Matthieu Constant's avatar
Matthieu Constant committed
39
use_lmgc = False
Matthieu Constant's avatar
Matthieu Constant committed
40
41
42
43
44

def genInitialPosition(filename, N, r, lx, ly, rhop) :

    p = scontact.ParticleProblem(2)
    #Loading of the mesh.msh file specifying physical boundaries name
Matthieu Constant's avatar
Matthieu Constant committed
45
    p.load_msh_boundaries("mesh.msh", ["Top", "Bottom", "Lateral"])
Matthieu Constant's avatar
Matthieu Constant committed
46
    #Definition of the points where the grains are located
Matthieu Constant's avatar
Matthieu Constant committed
47
48
    Ra = 0.8*r
    Rb = 0.9*r
Matthieu Constant's avatar
Matthieu Constant committed
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
    radii = pre.granulo_Random(N, Ra, Rb)
    [nb_laid_particles, coors] = pre.depositInBox2D(radii,lx,ly)
    for i in range(nb_laid_particles):
       p.add_particle(coors[2*i:2*i+2], radii[i], radii[i]**2 * np.pi * rhop);
    p.write_vtk(filename,0,0)

#physical parameters
alpha = 0*np.pi/4.
g =  -9.81*np.cos(alpha)                                      # gravity
print(g)
rho0 = 1.117                                      # fluid density
rho1 = 785.92
nu0 = 1.57e-5                                   # kinematic viscosity
nu1 = 1.2e-3/rho1
tEnd = 50                                     # final time
r = 5e-4/2
N = 100000
lx = .136
ly = .09
rhop = 1059
#numerical parameters
dt = .001                                  # time step

shutil.copy("mesh.msh", outputdir +"/mesh.msh")
Matthieu Constant's avatar
Matthieu Constant committed
73
#genInitialPosition(outputdir, N, r, lx, ly, rhop)
Matthieu Constant's avatar
inject    
Matthieu Constant committed
74
75
76
77
78
79
80
81
if use_lmgc:
    friction=0.3
    lmgc90Interface.scontactTolmgc90(outputdir, 2, 0, friction)
    p = lmgc90Interface.ParticleProblem(2)
    p.write_vtk(outputdir,0,0)
else:
    p = scontact.ParticleProblem(2)
    p.read_vtk(outputdir,0)
Matthieu Constant's avatar
inject    
Matthieu Constant committed
82

Matthieu Constant's avatar
Matthieu Constant committed
83
84
85
86
87
88
89
90
91
fluid = fluid.FluidProblem(2,g,[nu0*rho0,nu1*rho1],[rho0,rho1])
fluid.load_msh("mesh.msh")

Q = 5e-5/60
w = 2e-3
l = 8e-3
L = 0.136
sc_v = Q/(l*w)

Matthieu Constant's avatar
Matthieu Constant committed
92
# number of iterations between output files
93
outf = 1
Matthieu Constant's avatar
Matthieu Constant committed
94
95
96
ii = 0
t = 0
def outerBndV(x) :
Matthieu Constant's avatar
Matthieu Constant committed
97
98
99
100
    v = np.ndarray((len(x[:,0])))
    v[:] = 0
    v[np.abs(x[:,0]-L/2)<l/2] = sc_v
    return v
Matthieu Constant's avatar
Matthieu Constant committed
101

102
103
fluid.set_wall_boundary("Lateral",pressure=None)
fluid.set_open_boundary("Top",pressure=0,porosity=1,concentration=1)
Matthieu Constant's avatar
Matthieu Constant committed
104
fluid.set_open_boundary("Bottom",velocity=[0,outerBndV],porosity=1,concentration=1)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
105
#fluid.set_strong_boundary("Injection",1,outerBndV)
Matthieu Constant's avatar
Matthieu Constant committed
106

107
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
Matthieu Constant's avatar
Matthieu Constant committed
108
109

#set initial_condition
Matthieu Constant's avatar
Matthieu Constant committed
110
111
c = np.zeros_like(fluid.coordinates()[:,1])
fluid.set_concentration_cg(c)
Matthieu Constant's avatar
Matthieu Constant committed
112
113
114
115
116
117

fluid.export_vtk(outputdir,0,0)

ii = 0
tic = time.time() 

118
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
Matthieu Constant's avatar
Matthieu Constant committed
119
120
121
122
while t < tEnd : 
#Adaptation of the mesh. Args are minimal mesh radius, maximal mesh radius and number of elements
    #Fluid solver
    fluid.implicit_euler(dt, newton_max_it=20)
Matthieu Constant's avatar
Matthieu Constant committed
123
124
    #if (ii%50==0):
     #  fluid.adapt_mesh(1e-2,1e-2/5,5000)
Matthieu Constant's avatar
Matthieu Constant committed
125
126
127
128
129
    forces = fluid.compute_node_force(dt)
    #Computation of the new velocities
    vn = p.velocity() + forces * dt / p.mass()
    vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
    #number of sub time step
Matthieu Constant's avatar
Matthieu Constant committed
130
    nsub = max(5, int(np.ceil((vmax * dt * 4)/min(p.r()))))
Matthieu Constant's avatar
Matthieu Constant committed
131
132
133
134
135
    print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
    #Contact solver
    for i in range(nsub) :
        p.iterate(dt/nsub, forces)  
    t += dt
136
    fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
Matthieu Constant's avatar
Matthieu Constant committed
137
138
139
140
141
142
143
    #Output files writting
    if ii %outf == 0 :
        ioutput = int(ii/outf) + 1
        p.write_vtk(outputdir, ioutput, t)
        fluid.export_vtk(outputdir, t, ioutput)
    ii += 1
    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))