Skip to content
GitLab
Projects
Groups
Snippets
/
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
fluidparticles
MigFlow
Commits
8cd8c59f
Commit
8cd8c59f
authored
Oct 30, 2018
by
Matthieu Constant
Browse files
couette comments
parent
a7519b24
Pipeline
#4500
failed with stage
in 14 seconds
Changes
2
Pipelines
1
Hide whitespace changes
Inline
Side-by-side
testcases/couette-2d/betonlmgc.py
View file @
8cd8c59f
...
...
@@ -21,15 +21,15 @@
#!/usr/bin/env python
'''
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.
'''
# TESTCASE DESCRIPTION
#
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.
#variable to use lmgc or not
use_lmgc
=
1
from
migflow
import
fluid
as
fluid
from
migflow
import
scontact
2
from
migflow
import
fluid
from
migflow
import
scontact
if
use_lmgc
:
from
migflow
import
lmgc2Interface
import
numpy
as
np
...
...
@@ -37,155 +37,147 @@ import os
import
time
import
shutil
'''
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
'''
def
genInitialPosition
(
filename
,
r
,
rout
,
rhop
)
:
#Grains structure builder
p
=
scontact2
.
ParticleProblem
()
#Load mesh.msh file specifying physical boundaries names
"""Create a container to initialize particles and write particles in an initial output file.
Keyword arguments:
filename -- directory to write the output file
r -- max particle radius
rout -- outer radius of the drum
rhop -- particle density
"""
# Grains structure builder
p
=
scontact
.
ParticleProblem
(
2
)
# Load mesh.msh file specifying physical boundaries names
p
.
load_msh_boundaries
(
"mesh.msh"
,
[
"Outer"
,
"Inner"
])
'''
First layer of grains
'''
#Definition of the points where the grains are located
x
=
np
.
arange
(
rout
,
-
rout
,
2
*
-
r
)
y
=
np
.
arange
(
-
rout
,
-
2
*
rout
/
3.
,
2
*
r
)
x
,
y
=
np
.
meshgrid
(
x
,
y
)
R2
=
x
**
2
+
y
**
2
#
c
ondition to be inside the outer boundary
keep
=
R2
<
(
rout
-
r
)
**
2
x
=
x
[
keep
]
y
=
y
[
keep
]
#Add a grain at each centre position
#
First layer of grains
#
Definition of the points where the grains are located
x
=
np
.
arange
(
rout
,
-
rout
,
2
*
-
r
)
y
=
np
.
arange
(
-
rout
,
-
2
*
rout
/
3.
,
2
*
r
)
x
,
y
=
np
.
meshgrid
(
x
,
y
)
R2
=
x
**
2
+
y
**
2
#
C
ondition to be inside the outer boundary
keep
=
R2
<
(
rout
-
r
)
**
2
x
=
x
[
keep
]
y
=
y
[
keep
]
#
Add a grain at each centre position
for
i
in
range
(
x
.
shape
[
0
])
:
p
.
add_particle
((
x
[
i
],
y
[
i
]),
r
,
r
**
2
*
np
.
pi
*
rhop
);
'''
Second layer of grains
'''
#Definition of the points where the grains are located
x
=
np
.
arange
(
rout
,
-
rout
,
-
r
)
y
=
np
.
arange
(
-
2
*
rout
/
3.
+
r
,
-
rout
/
2.
,
r
)
x
,
y
=
np
.
meshgrid
(
x
,
y
)
#
c
ondition to be inside the outer boundary
R2
=
x
**
2
+
y
**
2
keep
=
R2
<
(
rout
-
r
)
**
2
x
=
x
[
keep
]
y
=
y
[
keep
]
#Add a grain at each centre position
#
Second layer of grains
#
Definition of the points where the grains are located
x
=
np
.
arange
(
rout
,
-
rout
,
-
r
)
y
=
np
.
arange
(
-
2
*
rout
/
3.
+
r
,
-
rout
/
2.
,
r
)
x
,
y
=
np
.
meshgrid
(
x
,
y
)
#
C
ondition to be inside the outer boundary
R2
=
x
**
2
+
y
**
2
keep
=
R2
<
(
rout
-
r
)
**
2
x
=
x
[
keep
]
y
=
y
[
keep
]
#
Add a grain at each centre position
for
i
in
range
(
x
.
shape
[
0
])
:
p
.
add_particle
((
x
[
i
],
y
[
i
]),
r
/
2.
,
r
**
2
/
4.
*
np
.
pi
*
rhop
);
'''
Third layer of grains
'''
#Definition of the points where the grains are located
x
=
np
.
arange
(
rout
,
-
rout
,
2
*
-
r
/
3
)
y
=
np
.
arange
(
-
rout
/
2.
+
r
,
-
5
*
rout
/
12.
,
2
*
r
/
3
)
x
,
y
=
np
.
meshgrid
(
x
,
y
)
R2
=
x
**
2
+
y
**
2
#
c
ondition to be inside the outer boundary
keep
=
R2
<
(
rout
-
r
/
3
)
**
2
x
=
x
[
keep
]
y
=
y
[
keep
]
#Add a grain at each centre position
#
Third layer of grains
#
Definition of the points where the grains are located
x
=
np
.
arange
(
rout
,
-
rout
,
2
*
-
r
/
3
)
y
=
np
.
arange
(
-
rout
/
2.
+
r
,
-
5
*
rout
/
12.
,
2
*
r
/
3
)
x
,
y
=
np
.
meshgrid
(
x
,
y
)
R2
=
x
**
2
+
y
**
2
#
C
ondition to be inside the outer boundary
keep
=
R2
<
(
rout
-
r
/
3
)
**
2
x
=
x
[
keep
]
y
=
y
[
keep
]
#
Add a grain at each centre position
for
i
in
range
(
x
.
shape
[
0
])
:
p
.
add_particle
((
x
[
i
],
y
[
i
]),
r
/
3
,
r
**
2
/
9
*
np
.
pi
*
rhop
);
p
.
add_particle
((
x
[
i
],
y
[
i
]),
r
/
3
,
r
**
2
/
9
*
np
.
pi
*
rhop
);
p
.
write_vtk
(
filename
,
0
,
0
)
#Define output directory
#
Define output directory
outputdir
=
"output"
if
not
os
.
path
.
isdir
(
outputdir
)
:
os
.
makedirs
(
outputdir
)
#
p
hysical parameters
g
=
-
9.81
# gravity
rhop
=
1600
# grains density
r
=
397e-6
/
2
# grains radius
rho
=
1e3
# fluid density
nu
=
1e-4
# kinematic viscosity
#
P
hysical parameters
g
=
-
9.81
# gravity
rhop
=
1600
# grains density
r
=
397e-6
/
2
# grains radius
rho
=
1e3
# fluid density
nu
=
1e-4
# kinematic viscosity
#
n
umerical parameters
outf
=
1
# number of iterations between output files
dt
=
5e-3
# time step
tEnd
=
100
# final time
#
N
umerical parameters
outf
=
1
# number of iterations between output files
dt
=
5e-3
# time step
tEnd
=
100
# final time
#Geometrical parametres
rout
=
0.0254
# outer boundary
#
Geometrical parametres
rout
=
0.0254
# outer boundary
#Initialise grains
#
# PARTICLE PROBLEM
#
# Initialise grains
genInitialPosition
(
outputdir
,
r
,
rout
,
rhop
)
if
use_lmgc
:
friction
=
0.1
# friction coefficient
lmgc2Interface
.
scontactToLmgc90
(
outputdir
,
0
,
friction
)
p
=
lmgc2Interface
.
LmgcInterface
()
p
=
lmgc2Interface
.
LmgcInterface
()
else
:
p
=
scontact
2
.
ParticleProblem
()
p
=
scontact
.
ParticleProblem
(
2
)
p
.
read_vtk
(
outputdir
,
0
)
#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
])
# Initial time and iteration
t
=
0
ii
=
0
#
# FLUID PROBLEM
#
fluid
=
fluid
.
FluidProblem
(
2
,
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
:
-
4
*
x
[:,
1
])
fluid
.
set_strong_boundary
(
"Outer"
,
1
,
lambda
x
:
4
*
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"
)
#Set location of the grains in the mesh and compute the porosity in each computation cell
fluid
.
set_particles
(
p
.
mass
(),
p
.
volume
(),
p
.
position
(),
p
.
velocity
())
fluid
.
export_vtk
(
outputdir
,
0
,
0
)
tic
=
time
.
time
()
#Computation loop
tic
=
time
.
time
()
fluid
.
set_particles
(
p
.
mass
(),
p
.
volume
(),
p
.
position
(),
p
.
velocity
())
#
# COMPUTATION LOOP
#
while
t
<
tEnd
:
#Fluid solver
#
Fluid solver
fluid
.
implicit_euler
(
dt
)
#Computation of the new velocities
forces
=
fluid
.
compute_node_force
(
dt
)
vn
=
p
.
velocity
()
+
forces
*
dt
/
p
.
mass
()
vmax
=
np
.
max
(
np
.
hypot
(
vn
[:,
0
],
vn
[:,
1
]))
#
n
umber of contact sub time step
nsub
=
max
(
1
,
int
(
np
.
ceil
((
vmax
*
dt
*
4
)
/
min
(
p
.
r
()))))
#
Computation of the new velocities
forces
=
fluid
.
compute_node_force
(
dt
)
vn
=
p
.
velocity
()
+
forces
*
dt
/
p
.
mass
()
vmax
=
np
.
max
(
np
.
hypot
(
vn
[:,
0
],
vn
[:,
1
]))
#
N
umber of contact 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
#
NLGS iterations
for
i
in
range
(
nsub
)
:
tol
=
1e-6
#numerical tolerance for grains intersection
tol
=
1e-6
#numerical tolerance for grains intersection
p
.
iterate
(
dt
/
nsub
,
forces
,
tol
)
#Localisation of the grains in the fluid
#
Localisation of the grains in the fluid
fluid
.
set_particles
(
p
.
mass
(),
p
.
volume
(),
p
.
position
(),
p
.
velocity
())
t
+=
dt
t
+=
dt
#Output files writting
#
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
))
ii
+=
1
print
(
"%i : %.2g/%.2g (cpu %.6g)"
%
(
ii
,
t
,
tEnd
,
time
.
time
()
-
tic
))
\ No newline at end of file
testcases/couette-2d/melangeur.py
View file @
8cd8c59f
...
...
@@ -20,13 +20,14 @@
# see <http://www.gnu.org/licenses/>.
#!/usr/bin/env python
'''
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.
'''
from
migflow
import
fluid
as
fluid
from
migflow
import
scontact2
# TESTCASE DESCRIPTION
# 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.
from
migflow
import
fluid
from
migflow
import
scontact
import
numpy
as
np
import
os
...
...
@@ -34,129 +35,120 @@ import time
import
shutil
import
random
#Define output directory
#
Define output directory
outputdir
=
"output"
if
not
os
.
path
.
isdir
(
outputdir
)
:
os
.
makedirs
(
outputdir
)
'''
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
'''
def
genInitialPosition
(
filename
,
r
,
rout
,
rin
,
rhop
)
:
#Grains structure builder
p
=
scontact2
.
ParticleProblem
()
#Load mesh.msh file specifying physical boundaries names
"""Create a container to initialize particles and write particles in an initial output file.
Keyword arguments:
filename -- directory to write the output file
r -- max particle radius
rout -- outer radius of the drum
rhop -- particle density
"""
# Grains structure builder
p
=
scontact
.
ParticleProblem
(
2
)
# Load mesh.msh file specifying physical boundaries names
p
.
load_msh_boundaries
(
"mesh.msh"
,
[
"Outer"
,
"Inner"
])
#Definition of the points where the grains are located
x
=
np
.
arange
(
rout
,
-
rout
,
2
*
-
r
)
x
,
y
=
np
.
meshgrid
(
x
,
x
)
R2
=
x
**
2
+
y
**
2
#
c
ondition to be inside the outer boundary
keep
=
R2
<
(
rout
-
r
)
**
2
x
=
x
[
keep
]
y
=
y
[
keep
]
R2
=
x
**
2
+
y
**
2
#
c
ondition to be outside the inner boundary
keep
=
R2
>
(
rin
+
r
)
**
2
x
=
x
[
keep
]
y
=
y
[
keep
]
#Add a grain at each centre position
#
Definition of the points where the grains are located
x
=
np
.
arange
(
rout
,
-
rout
,
2
*
-
r
)
x
,
y
=
np
.
meshgrid
(
x
,
x
)
R2
=
x
**
2
+
y
**
2
#
C
ondition to be inside the outer boundary
keep
=
R2
<
(
rout
-
r
)
**
2
x
=
x
[
keep
]
y
=
y
[
keep
]
R2
=
x
**
2
+
y
**
2
#
C
ondition to be outside the inner boundary
keep
=
R2
>
(
rin
+
r
)
**
2
x
=
x
[
keep
]
y
=
y
[
keep
]
#
Add a grain at each centre position
for
i
in
range
(
x
.
shape
[
0
])
:
if
y
[
i
]
<
0
:
rhop1
=
random
.
choice
([
rhop
*
.
9
,
1.1
*
rhop
,
rhop
])
p
.
add_particle
((
x
[
i
],
y
[
i
]),
r
,
r
**
2
*
np
.
pi
*
rhop1
)
p
.
write_vtk
(
filename
,
0
,
0
)
#physical parameters
g
=
0
# gravity
rho
=
1.253e3
# fluid density
rhop
=
1000
# grains density
nu
=
1e-4
# kinematic viscosity
#numerical parameters
dt
=
1e-2
# time step
tEnd
=
50
# final time
#geometry parameters
outf
=
5
# number of iterations between output files
rout
=
0.0254
# outer radius
rin
=
0.0064
# inner radius
r
=
397e-6
/
2
# grains radius
#Object particles creation
# Physical parameters
g
=
0
# gravity
rho
=
1.253e3
# fluid density
rhop
=
1000
# grains density
nu
=
1e-4
# kinematic viscosity
# Numerical parameters
dt
=
1e-2
# time step
tEnd
=
50
# final time
# Geometry parameters
outf
=
5
# number of iterations between output files
rout
=
0.0254
# outer radius
rin
=
0.0064
# inner radius
r
=
397e-6
/
2
# grains radius
#
# PARTICLE PROBLEM
#
# Object particles creation
genInitialPosition
(
outputdir
,
r
,
rout
,
rin
,
rhop
)
p
=
scontact
2
.
ParticleProblem
()
p
=
scontact
.
ParticleProblem
(
2
)
p
.
read_vtk
(
outputdir
,
0
)
#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
# Initial time and iteration
t
=
0
ii
=
0
#
# FLUID PROBLEM
#
fluid
=
fluid
.
FluidProblem
(
2
,
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"
)
#Set location of the grains in the mesh and compute the porosity in each computation cell
# Set location of the grains in the mesh and compute the porosity in each computation cell
fluid
.
set_particles
(
p
.
mass
(),
p
.
volume
(),
p
.
position
(),
p
.
velocity
())
fluid
.
export_vtk
(
outputdir
,
0
,
0
)
tic
=
time
.
time
()
#Computation loop
fluid
.
set_particles
(
p
.
mass
(),
p
.
volume
(),
p
.
position
(),
p
.
velocity
())
#
# COMPUTATION LOOP
#
while
t
<
tEnd
:
#Fluid solver
#
Fluid solver
fluid
.
implicit_euler
(
dt
)
#Computation of the new velocities
forces
=
fluid
.
compute_node_force
(
dt
)
vn
=
p
.
velocity
()
+
forces
*
dt
/
p
.
mass
()
vmax
=
np
.
max
(
np
.
hypot
(
vn
[:,
0
],
vn
[:,
1
]))
#
n
umber of contact sub time step
nsub
=
max
(
1
,
int
(
np
.
ceil
((
vmax
*
dt
*
4
)
/
min
(
p
.
r
()))))
#
Computation of the new velocities
forces
=
fluid
.
compute_node_force
(
dt
)
vn
=
p
.
velocity
()
+
forces
*
dt
/
p
.
mass
()
vmax
=
np
.
max
(
np
.
hypot
(
vn
[:,
0
],
vn
[:,
1
]))
#
N
umber of contact 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
#
NLGS iterations
for
i
in
range
(
nsub
)
:
tol
=
1e-6
#
n
umerical tolerance for grains intersection
tol
=
1e-6
#
N
umerical tolerance for grains intersection
p
.
iterate
(
dt
/
nsub
,
forces
,
tol
)
#Localisation of the grains in the fluid
#
Localisation of the grains in the fluid
fluid
.
set_particles
(
p
.
mass
(),
p
.
volume
(),
p
.
position
(),
p
.
velocity
())
t
+=
dt
t
+=
dt
#Output files writting
#
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
))
ii
+=
1
print
(
"%i : %.2g/%.2g (cpu %.6g)"
%
(
ii
,
t
,
tEnd
,
time
.
time
()
-
tic
))
\ No newline at end of file
Write
Preview
Supports
Markdown
0%
Try again
or
attach a new file
.
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Cancel
Please
register
or
sign in
to comment