dep.py 6.17 KB
Newer Older
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
1
# MigFlow - Copyright (C) <2010-2018>
2
3
# <Universite catholique de Louvain (UCL), Belgium
#  Universite de Montpellier, France>
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
4
# 	
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
5
# List of the contributors to the development of MigFlow: see AUTHORS file.
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
6
7
# Description and complete License: see LICENSE file.
# 	
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
8
# This program (MigFlow) is free software: 
9
# you can redistribute it and/or modify it under the terms of the GNU Lesser General 
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
10
11
12
13
14
15
# 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
16
# GNU Lesser General Public License for more details.
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
17
# 
18
19
# 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, 
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
20
21
# see <http://www.gnu.org/licenses/>.

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
22
#!/usr/bin/env python
Matthieu Constant's avatar
Matthieu Constant committed
23
24
25
26
27
28

'''
Mixing of grains in a rotating drum.
This example shows how to set boundary conditions as a function of some parametres.
A boolean parametre give the possibility to use lmgc90 to solve contacts or not.
'''
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
29
30
from migflow import fluid as fluid
from migflow import scontact2
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
31
32
33
34
35
36
37

import numpy as np
import os
import time
import shutil
import random

Matthieu Constant's avatar
Matthieu Constant committed
38
39
40
41
42
43
44
45
46
47
48
'''
genInitialPosition is a function that sets all the grains centre positions and create the grains objects to add in the computing structure
Args:    -filename    : name of the output file
         -r           : max radius of the grains
         -H           : domain height
         -ly          : grains area height
         -lx          : grains area width
         -rhop        : grains density        
'''
def genInitialPosition(filename, r, H, ly, lx, rhop) :
    #Grains structure builder
Matthieu Constant's avatar
Matthieu Constant committed
49
    p     = scontact2.ParticleProblem()
Matthieu Constant's avatar
Matthieu Constant committed
50
    #Load mesh.msh file specifying physical boundaries names
51
    p.load_msh_boundaries("mesh.msh", ["Top", "Lateral","Bottom"])
Matthieu Constant's avatar
Matthieu Constant committed
52

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
53
    #Definition of the points where the grains are located
Matthieu Constant's avatar
Matthieu Constant committed
54
55
56
57
58
    x     = np.arange(-lx/2+r, lx/2-r, 2*r)
    y     = np.arange(H/2-r, H/2-ly+r, -2*r)
    x, y  = np.meshgrid(x, y)
    x     = x.flat
    y     = y.flat
Matthieu Constant's avatar
Matthieu Constant committed
59
    #Add a grain at each centre position
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
60
    for i in range(len(x)) :
Matthieu Constant's avatar
Matthieu Constant committed
61
        p.add_particle((x[i], y[i]), r, r**2 * np.pi * rhop)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
62
63
    p.write_vtk(filename,0,0)

Matthieu Constant's avatar
Matthieu Constant committed
64
#Define output directory
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
65
66
67
68
69
70
71
72
73
outputdir = "output"
if not os.path.isdir(outputdir) :
    os.makedirs(outputdir)

t = 0
ii = 0


#physical parameters
Matthieu Constant's avatar
Matthieu Constant committed
74
75
76
77
78
g         = -9.81                                     # gravity
r         =  1e-3                                     # grains radius
rhop      =  1500                                     # grains density
rho       =  1000                                     # fluid density
nu        =  1e-6                                     # kinematic viscosity
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
79
80

#numerical parameters
Matthieu Constant's avatar
Matthieu Constant committed
81
outf      =  5                                        # number of iterations between output files
Matthieu Constant's avatar
Matthieu Constant committed
82
83
dt        =  1e-2                                     # time step
tEnd      =  100                                      # final time
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
84

Matthieu Constant's avatar
Matthieu Constant committed
85
86
87
88
#geometrical parameters
ly        =  5e-2                                     # grains area height
lx        =  4e-1                                     # grains area widht
H         =  1.2                                      # domain height
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
89

Matthieu Constant's avatar
Matthieu Constant committed
90
91
92
93

#Initialise grains
genInitialPosition(outputdir, r, H, ly, lx, rhop)
p         =  scontact2.ParticleProblem()
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
94
95
96
97
98
p.read_vtk(outputdir,0)

print("r = %g, m = %g\n" %  (p.r()[0], p.mass()[0]))
print("RHOP = %g" % rhop)

Matthieu Constant's avatar
Matthieu Constant committed
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
#Initial time and iteration
t         =  0
ii        =  0

'''
fluid_problem is the function builder of the fluid class
Args:     -g             : gravity
          -viscosity     : list of fluid phases viscosities
          -density       : list of fluid phases densities
'''
fluid = fluid.fluid_problem(g,[nu*rho],[rho])
#Set the mesh geometry for the fluid computation
fluid.load_msh("mesh.msh")
'''
set_strong_boundary is the function setting the strong boundaries (=constraints) for the fluid problem
Args:     -Tag           : physical tag (set in the mesh.geo file) of the boundary on which you want to add a constraint
          -Field         : O: x-velocity; 1: y-velocity; 2: pressure
          -Value         : value assigned to the specified field on the specified boundary
'''
fluid.set_strong_boundary("Top",2,0)
fluid.set_strong_boundary("Top",1,0)
fluid.set_strong_boundary("Bottom",1,0)
fluid.set_strong_boundary("Lateral",0,0)
'''
set_weak_boundary is the function setting the weak boundaries (=normal fluxes) for the fluid problem
Args:     -Tag           : physical tag (set in the mesh.geo file) of the boundary on which you want to add a constraint
          -Bnd Type      : boundaries can be of type: wall, inflow, outflow
'''
127
128
129
fluid.set_weak_boundary("Bottom","Wall")
fluid.set_weak_boundary("Lateral","Wall")
fluid.set_weak_boundary("Top","Wall")
Matthieu Constant's avatar
Matthieu Constant committed
130
#Set location of the grains in the mesh and compute the porosity in each computation cell
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
131
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
132
fluid.export_vtk(outputdir,0,0)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
133
134


Matthieu Constant's avatar
Matthieu Constant committed
135
tic = time.time()
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
136
137
138
139
140
#Computation loop
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
while t < tEnd :
    #Fluid solver
    fluid.implicit_euler(dt)
Matthieu Constant's avatar
Matthieu Constant committed
141

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
142
    #Computation of the new velocities
Matthieu Constant's avatar
Matthieu Constant committed
143
144
145
146
147
    forces      = fluid.compute_node_force(dt)
    vn          = p.velocity() + forces * dt / p.mass()
    vmax        = np.max(np.hypot(vn[:, 0], vn[:, 1]))
    #number of contact sub time step
    nsub        = max(1, int(np.ceil((vmax * dt * 4)/min(p.r()))))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
148
149
150
    print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
    #Contact solver
    for i in range(nsub) :
Matthieu Constant's avatar
Matthieu Constant committed
151
        tol     = 1e-6                                #numerical tolerance for grains intersection
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
152
        p.iterate(dt/nsub, forces, tol)
Matthieu Constant's avatar
Matthieu Constant committed
153

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
154
155
    #Localisation of the grains in the fluid
    fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
Matthieu Constant's avatar
Matthieu Constant committed
156
157
    t          += dt

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
158
159
160
161
162
    #Output files writting
    if ii %outf == 0 :
        ioutput = int(ii/outf) + 1
        p.write_vtk(outputdir, ioutput, t)
        fluid.export_vtk(outputdir, t, ioutput)
Matthieu Constant's avatar
Matthieu Constant committed
163
    ii         += 1
Matthieu Constant's avatar
Matthieu Constant committed
164
    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))