Commit 177fee02 authored by Matthieu Constant's avatar Matthieu Constant
Browse files

melangeur with comments

parent 29b60095
Pipeline #4477 passed with stage
in 1 minute and 10 seconds
......@@ -20,6 +20,11 @@
# 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
......@@ -29,100 +34,129 @@ import time
import shutil
import random
#Define output directory
outputdir = "output"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
#grains creation and placement + physical boundaries definition
'''r is the radius of the grains
rout is the outer radius of the geometry
rin is the inner radius of the geometry
rhop is the particles density'''
'''
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) :
p = scontact2.ParticleProblem()
#Loading of the mesh.msh file specifying physical boundaries name
#Grains structure builder
p = scontact2.ParticleProblem()
#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
x = np.arange(rout, -rout, 2 * -r)
x, y = np.meshgrid(x, x)
R2 = x**2 + y**2
#condition to be inside the outer boundary
keep = R2 < (rout - r)**2
x = x[keep]
y = y[keep]
R2 = x**2 + y**2
keep = R2 < (rout - r)**2
x = x[keep]
y = y[keep]
R2 = x**2 + y**2
#condition to be outside the inner boundary
keep = R2 > (rin + r)**2
x = x[keep]
y = y[keep]
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])
#Addition of an particle object at each point
p.add_particle((x[i], y[i]), r, r**2 * np.pi * rhop1);
p.write_vtk(filename,0,0)
t = 0
ii = 0
#physical parameters
g = 0 # gravity
rho = 1.253e3 # fluid density
rhop = 1000 # grains density
nu = 1e-4 # kinematic viscosity
tEnd = 50 # final time
g = 0 # gravity
rho = 1.253e3 # fluid density
rhop = 1000 # grains density
nu = 1e-4 # kinematic viscosity
#numerical parameters
h = 0.002 # mesh size
dt = 1e-2 # time step
epsilon = 1e-5 # stabilisation parametre
dt = 1e-2 # time step
tEnd = 50 # final time
#geometry parameters
rout = 0.0254 # outer radius
rin = 0.0064 # inner radius
r = 397e-6/2 # grains radius
shutil.copy("mesh.msh", outputdir +"/mesh.msh")
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
genInitialPosition(outputdir, r, rout, rin, rhop)
p = scontact2.ParticleProblem()
p = scontact2.ParticleProblem()
p.read_vtk(outputdir,0)
outf = 5 # number of iterations between output files
#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")
#Object fluid creation + Boundary condition of the fluid (field 0 is horizontal velocity; field 1 is vertical velocity; field 2 is pressure)
#Format: strong_boundaries = [(Boundary tag, Fluid field, Value)
strong_boundaries = [("Outer",0,lambda x : -x[:, 1]),("Outer",1,lambda x : x[:, 0]),("Inner",0,0.),("Inner",1,0.),("PtFix",2,0.)]
fluid = fluid.fluid_problem("mesh.msh",g,nu*rho,rho,epsilon,strong_boundaries)
#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(outputdir,0,0)
fluid.export_vtk(outputdir,0,0)
tic = time.clock()
fluid.export_vtk(outputdir, 0, 0s)
tic = time.time()
#Computation loop
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
while t < tEnd :
#Fluid solver
fluid.implicit_euler(dt)
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()))))
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()))))
print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
#Contact solver
for i in range(nsub) :
tol = 1e-6
tol = 1e-6 #numerical tolerance for grains intersection
p.iterate(dt/nsub, forces, tol)
#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
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
p.write_vtk(outputdir, ioutput, t)
fluid.export_vtk(outputdir, t, ioutput)
ii += 1
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.clock() - tic))
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment