melangeur.py 3.73 KB
Newer Older
1
2
3
4
5
6
7
8
#!/usr/bin/env python
from marblesbag import fluid as fluid
from marblesbag import scontact2

import numpy as np
import os
import time
import shutil
Matthieu Constant's avatar
Matthieu Constant committed
9
import random
10
11
12
13
14

outputdir = "output"
if not os.path.isdir(outputdir) :
    os.makedirs(outputdir)

Matthieu Constant's avatar
Matthieu Constant committed
15
16
17
18
19
#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'''
Matthieu Constant's avatar
Matthieu Constant committed
20
def genInitialPosition(filename, r, rout, rin, rhop) :
21
    p = scontact2.ParticleProblem()
Matthieu Constant's avatar
Matthieu Constant committed
22
    #Loading of the mesh.msh file specifying physical boundaries name
23
    p.load_msh_boundaries("mesh.msh", ["Outer", "Inner"])
Matthieu Constant's avatar
Matthieu Constant committed
24
    #Definition of the points where the grains are located
25
26
27
    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
28
    #condition to be inside the outer boundary
29
30
31
    keep = R2 < (rout - r)**2
    x = x[keep]
    y = y[keep]
Matthieu Constant's avatar
Matthieu Constant committed
32
    R2 = x**2 + y**2
Matthieu Constant's avatar
Matthieu Constant committed
33
    #condition to be outside the inner boundary
Matthieu Constant's avatar
Matthieu Constant committed
34
35
36
    keep = R2 > (rin + r)**2
    x = x[keep]
    y = y[keep]
37
    for i in range(x.shape[0]) :
Matthieu Constant's avatar
Matthieu Constant committed
38
39
40
41
        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);
Matthieu Constant's avatar
Matthieu Constant committed
42
    p.write_vtk(filename,0,0)
43
44
45
46
47
48


t = 0
ii = 0

#physical parameters
Matthieu Constant's avatar
Matthieu Constant committed
49
50
51
52
53
g =  0                                  # gravity
rho = 1.253e3                           # fluid density
rhop = 1000                             # grains density
nu = 1e-4                               # kinematic viscosity
tEnd = 50                               # final time
54
55

#numerical parameters
Matthieu Constant's avatar
Matthieu Constant committed
56
57
58
h = 0.002                               # mesh size
dt = 1e-2                               # time step
epsilon = 1e-5                          # stabilisation parametre
59

Matthieu Constant's avatar
Matthieu Constant committed
60
61
62
63
#geometry parameters
rout = 0.0254                           # outer radius
rin = 0.0064                            # inner radius
r = 397e-6/2                            # grains radius
64
shutil.copy("mesh.msh", outputdir +"/mesh.msh")
Matthieu Constant's avatar
Matthieu Constant committed
65
66

#Object particles creation
Matthieu Constant's avatar
Matthieu Constant committed
67
68
69
genInitialPosition(outputdir, r, rout, rin, rhop)
p = scontact2.ParticleProblem()
p.read_vtk(outputdir,0)
70

Matthieu Constant's avatar
Matthieu Constant committed
71
outf = 5                                # number of iterations between output files
72

Matthieu Constant's avatar
Matthieu Constant committed
73
74
#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)
75
strong_boundaries = [("Outer",0,lambda x : -x[:, 1]),("Outer",1,lambda x : x[:, 0]),("Inner",0,0.),("Inner",1,0.),("PtFix",2,0.)]
Matthieu Constant's avatar
Matthieu Constant committed
76
fluid = fluid.fluid_problem("mesh.msh",g,nu*rho,rho,epsilon,strong_boundaries)
77
78
79
80
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.export(outputdir,0,0)

tic = time.clock()
Matthieu Constant's avatar
Matthieu Constant committed
81
#Computation loop
Matthieu Constant's avatar
Matthieu Constant committed
82
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
Matthieu Constant's avatar
Matthieu Constant committed
83
84
while t < tEnd :
    #Fluid solver
Matthieu Constant's avatar
Matthieu Constant committed
85
    fluid.implicit_euler(dt)
86
    forces = fluid.compute_node_force(dt)
Matthieu Constant's avatar
Matthieu Constant committed
87
    #Computation of the new velocities
88
89
    vn = p.velocity() + forces * dt / p.mass()
    vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
Matthieu Constant's avatar
Matthieu Constant committed
90
    #number of sub time step
Matthieu Constant's avatar
Matthieu Constant committed
91
    nsub = max(1, int(np.ceil((vmax * dt * 4)/min(p.r()))))
92
    print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
Matthieu Constant's avatar
Matthieu Constant committed
93
    #Contact solver
94
    for i in range(nsub) :
Matthieu Constant's avatar
Matthieu Constant committed
95
96
97
        tol = 1e-6
        p.iterate(dt/nsub, forces, tol)
    #Localisation of the grains in the fluid
98
99
    fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
    t += dt
Matthieu Constant's avatar
Matthieu Constant committed
100
    #Output files writting
101
102
103
    if ii %outf == 0 :
        ioutput = int(ii/outf) + 1
        p.write_vtk(outputdir, ioutput, t)
Matthieu Constant's avatar
Matthieu Constant committed
104
        fluid.export_vtk(outputdir, t, ioutput)
105
106
    ii += 1
    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.clock() - tic))