melangeur.py 6.23 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/>.

22
#!/usr/bin/env python
Matthieu Constant's avatar
Matthieu Constant committed
23
24
25
26
27
'''
Grains in circular shear flow.
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
28
29
from migflow import fluid as fluid
from migflow import scontact2
30
31
32
33
34

import numpy as np
import os
import time
import shutil
Matthieu Constant's avatar
Matthieu Constant committed
35
import random
36

Matthieu Constant's avatar
Matthieu Constant committed
37
#Define output directory
38
39
40
41
outputdir = "output"
if not os.path.isdir(outputdir) :
    os.makedirs(outputdir)

Matthieu Constant's avatar
Matthieu Constant committed
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
         -rout        : outer radius of the drum
         -rhop        : grains density        
'''
Matthieu Constant's avatar
Matthieu Constant committed
49
def genInitialPosition(filename, r, rout, rin, rhop) :
Matthieu Constant's avatar
Matthieu Constant committed
50
51
52
    #Grains structure builder
    p     = scontact2.ParticleProblem()
    #Load mesh.msh file specifying physical boundaries names
53
    p.load_msh_boundaries("mesh.msh", ["Outer", "Inner"])
Matthieu Constant's avatar
Matthieu Constant committed
54
    
Matthieu Constant's avatar
Matthieu Constant committed
55
    #Definition of the points where the grains are located
Matthieu Constant's avatar
Matthieu Constant committed
56
57
58
    x     = np.arange(rout, -rout, 2 * -r)
    x, y  = np.meshgrid(x, x)
    R2    = x**2 + y**2
Matthieu Constant's avatar
Matthieu Constant committed
59
    #condition to be inside the outer boundary
Matthieu Constant's avatar
Matthieu Constant committed
60
61
62
63
    keep  = R2 < (rout - r)**2
    x     = x[keep]
    y     = y[keep]
    R2    = x**2 + y**2
Matthieu Constant's avatar
Matthieu Constant committed
64
    #condition to be outside the inner boundary
Matthieu Constant's avatar
Matthieu Constant committed
65
66
67
68
    keep  = R2 > (rin + r)**2
    x     = x[keep]
    y     = y[keep]
    #Add a grain at each centre position
69
    for i in range(x.shape[0]) :
Matthieu Constant's avatar
Matthieu Constant committed
70
71
        if y[i]<0:
            rhop1 = random.choice([rhop*.9,1.1*rhop,rhop])
Matthieu Constant's avatar
Matthieu Constant committed
72
            p.add_particle((x[i], y[i]), r, r**2 * np.pi * rhop1)
Matthieu Constant's avatar
Matthieu Constant committed
73
    p.write_vtk(filename,0,0)
74

Matthieu Constant's avatar
Matthieu Constant committed
75
 
76
#physical parameters
Matthieu Constant's avatar
Matthieu Constant committed
77
78
79
80
g         =  0                                        # gravity
rho       =  1.253e3                                  # fluid density
rhop      =  1000                                     # grains density
nu        =  1e-4                                     # kinematic viscosity
81
82

#numerical parameters
Matthieu Constant's avatar
Matthieu Constant committed
83
84
dt        =  1e-2                                     # time step
tEnd      =  50                                       # final time
85

Matthieu Constant's avatar
Matthieu Constant committed
86
#geometry parameters
Matthieu Constant's avatar
Matthieu Constant committed
87
88
89
90
outf      =  5                                        # number of iterations between output files
rout      =  0.0254                                   # outer radius
rin       =  0.0064                                   # inner radius
r         =  397e-6/2                                 # grains radius
Matthieu Constant's avatar
Matthieu Constant committed
91
92

#Object particles creation
Matthieu Constant's avatar
Matthieu Constant committed
93
genInitialPosition(outputdir, r, rout, rin, rhop)
Matthieu Constant's avatar
Matthieu Constant committed
94
p         =  scontact2.ParticleProblem()
Matthieu Constant's avatar
Matthieu Constant committed
95
p.read_vtk(outputdir,0)
96

Matthieu Constant's avatar
Matthieu Constant committed
97
98
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
127
#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("Outer",0,lambda x : -x[:, 1])
fluid.set_strong_boundary("Outer",1,lambda x : x[:, 0])
fluid.set_strong_boundary("Inner",0,0)
fluid.set_strong_boundary("Inner",1,0)
fluid.set_strong_boundary("PtFix",2,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
'''
fluid.set_weak_boundary("Inner","Wall")
fluid.set_weak_boundary("Outer","Wall")
128

Matthieu Constant's avatar
Matthieu Constant committed
129
#Set location of the grains in the mesh and compute the porosity in each computation cell
130
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
Matthieu Constant's avatar
Matthieu Constant committed
131
fluid.export_vtk(outputdir,0,0)
132

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

Matthieu Constant's avatar
Matthieu Constant committed
140
    #Computation of the new velocities
Matthieu Constant's avatar
Matthieu Constant committed
141
142
143
144
145
    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()))))
146
    print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
Matthieu Constant's avatar
Matthieu Constant committed
147
    #Contact solver
148
    for i in range(nsub) :
Matthieu Constant's avatar
Matthieu Constant committed
149
        tol     = 1e-6                                #numerical tolerance for grains intersection
Matthieu Constant's avatar
Matthieu Constant committed
150
        p.iterate(dt/nsub, forces, tol)
Matthieu Constant's avatar
Matthieu Constant committed
151

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

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