inject2fGrains.py 4.34 KB
Newer Older
Matthieu Constant's avatar
Matthieu Constant committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
# Marblesbag - Copyright (C) <2010-2018>
# <Universite catholique de Louvain (UCL), Belgium
#  Universite de Montpellier, France>
# 	
# List of the contributors to the development of Marblesbag: see AUTHORS file.
# Description and complete License: see LICENSE file.
# 	
# This program (Marblesbag) 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 marblesbag import fluid as fluid
from marblesbag import scontact2
Matthieu Constant's avatar
Matthieu Constant committed
25
26
from marblesbag import lmgc2Interface
from pylmgc90 import pre
Matthieu Constant's avatar
Matthieu Constant committed
27
28
29
30
31
32
33
import numpy as np
import os
import time
import shutil
import random


Matthieu Constant's avatar
Matthieu Constant committed
34
35


Matthieu Constant's avatar
Matthieu Constant committed
36
outputdir = "outputVid"
Matthieu Constant's avatar
Matthieu Constant committed
37
38
39
40
41
42
if not os.path.isdir(outputdir) :
    os.makedirs(outputdir)

t = 0
ii = 0

Matthieu Constant's avatar
Matthieu Constant committed
43
44
45
46
47
48
49
50
51
52
53
54
55
def genInitialPosition(filename, N, r, lx, ly, rhop) :

    p = scontact2.ParticleProblem()
    #Loading of the mesh.msh file specifying physical boundaries name
    p.load_msh_boundaries("mesh.msh", ["Top", "Bottom", "Lateral", "Injection"])
    #Definition of the points where the grains are located
    Ra = 0.9*r
    Rb = 1.1*r
    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)
Matthieu Constant's avatar
Matthieu Constant committed
56
57

#physical parameters
Matthieu Constant's avatar
Matthieu Constant committed
58
59
60
alpha = 0*np.pi/4.
g =  -9.81*np.cos(alpha)                                      # gravity
print(g)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
61
rho0 = 1.117                                      # fluid density
Matthieu Constant's avatar
Matthieu Constant committed
62
rho1 = 1000
Matthieu Constant's avatar
Matthieu Constant committed
63
64
nu1 = 1.57e-5                                   # kinematic viscosity
nu0 = 1e-6
Matthieu Constant's avatar
Matthieu Constant committed
65
66
tEnd = 10                                     # final time
r = 5e-4
Matthieu Constant's avatar
Matthieu Constant committed
67
N = 10000*3
Matthieu Constant's avatar
Matthieu Constant committed
68
lx = .222
Matthieu Constant's avatar
Matthieu Constant committed
69
ly = .14*2
Matthieu Constant's avatar
Matthieu Constant committed
70
rhop = 2500
Matthieu Constant's avatar
Matthieu Constant committed
71
#numerical parameters
Matthieu Constant's avatar
Matthieu Constant committed
72
dt = .0001                                  # time step
Matthieu Constant's avatar
Matthieu Constant committed
73
74

shutil.copy("mesh.msh", outputdir +"/mesh.msh")
Matthieu Constant's avatar
Matthieu Constant committed
75
76
genInitialPosition(outputdir, N, r, lx, ly, rhop)
friction=0.3
Matthieu Constant's avatar
Matthieu Constant committed
77
78
p = scontact2.ParticleProblem()
p.read_vtk(outputdir,0)
Matthieu Constant's avatar
Matthieu Constant committed
79
# number of iterations between output files
Matthieu Constant's avatar
Matthieu Constant committed
80
outf = 25
Matthieu Constant's avatar
Matthieu Constant committed
81
82
ii = 0
t = 0
Matthieu Constant's avatar
Matthieu Constant committed
83
def outerBndV(x) :
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
84
    return max(1*t,0.1)
Matthieu Constant's avatar
Matthieu Constant committed
85

86
fluid = fluid.fluid_problem(g,[nu0*rho0,nu1*rho1],[rho0,rho1])
87
fluid.load_msh("mesh.msh")
88
89
90
91
92
93

fluid.set_strong_boundary("Top",2,2,0)
fluid.set_strong_boundary("Injection",1,1,outerBndV)
fluid.set_strong_boundary("Injection",3,3,1)
fluid.set_strong_boundary("Injection",0,0,0)

Matthieu Constant's avatar
Matthieu Constant committed
94
95
96
97
fluid.set_weak_boundary("Bottom","Wall")
fluid.set_weak_boundary("Lateral","Wall")
fluid.set_weak_boundary("Top","Outflow")
fluid.set_weak_boundary("Injection","Inflow")
Matthieu Constant's avatar
Matthieu Constant committed
98

Matthieu Constant's avatar
Matthieu Constant committed
99
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
100

Matthieu Constant's avatar
Matthieu Constant committed
101
#set initial_condition
Matthieu Constant's avatar
Matthieu Constant committed
102
103
print(fluid.solution().shape, fluid.coordinates().shape)
fluid.solution()[fluid.coordinates()[:,1]>0.3,3] = 1
Matthieu Constant's avatar
Matthieu Constant committed
104
105
106
107
108

fluid.export_vtk(outputdir,0,0)

ii = 0
tic = time.clock()
Matthieu Constant's avatar
Matthieu Constant committed
109
110

fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
Matthieu Constant's avatar
Matthieu Constant committed
111
112
113
114
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)
Matthieu Constant's avatar
Matthieu Constant committed
115
    if (ii%65==0 and ii != 0):
Matthieu Constant's avatar
Matthieu Constant committed
116
       fluid.adapt_mesh(8e-2,3e-3,10000)
Matthieu Constant's avatar
Matthieu Constant committed
117
118
119
120
121
122
123
124
125
126
127
    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
    nsub = max(1, int(np.ceil((vmax * dt * 4)/min(p.r()))))
    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
Matthieu Constant's avatar
Matthieu Constant committed
128
    fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
Matthieu Constant's avatar
Matthieu Constant committed
129
130
131
132
133
134
135
    #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.clock() - tic))