depot.py 3.98 KB
Newer Older
Matthieu Constant's avatar
Matthieu Constant committed
1
2
3
4
5
6
7
8
9
#!/usr/bin/env python
from marblesbag import fluid3 as fluid
from marblesbag import scontact3

import numpy as np
import os
import time
import shutil

Matthieu Constant's avatar
Matthieu Constant committed
10
11
#grains creation and placement + physical boundaries definition
#The radii of the grains are defined in radius.csv and are computed to fit the distribution size given in the article of Andre et al. (2011) "Real-time analysis of the growth of granular media by an ultrasonic method: Application to the sedimentation of glass balls in water"
Matthieu Constant's avatar
Matthieu Constant committed
12
13
def genInitialPosition(filename, rhop) :
    p = scontact3.ParticleProblem()
Matthieu Constant's avatar
Matthieu Constant committed
14
    #Loading of the mesh.msh file specifying physical boundaries name
Matthieu Constant's avatar
Matthieu Constant committed
15
    p.load_msh_boundaries("mesh.msh", ["Top", "Box"])
Matthieu Constant's avatar
Matthieu Constant committed
16
    #Reading the file containing the radius
Matthieu Constant's avatar
Matthieu Constant committed
17
18
19
    myRadius = np.genfromtxt('radius.csv', delimiter=',')
    r2=np.amax(myRadius)*1e-6
    i = 0
Matthieu Constant's avatar
Matthieu Constant committed
20
    #Positioning the grains on a regular grid
Matthieu Constant's avatar
Matthieu Constant committed
21
22
23
24
25
    for y in np.arange(0, -.004, -2*r2):
        for z in np.arange(.01675-r2, -.01675, -2*r2):
            for x in np.arange(0.014-r2, -.014+r2, -2*r2):
                p.add_particle((x, y, z), myRadius[i%500000]*1e-6, 4./3.* (myRadius[i%500000]*1e-6)**3 * np.pi * rhop);            
                i += 1
Matthieu Constant's avatar
Matthieu Constant committed
26
    #Translation of the grains to the top of the box
Matthieu Constant's avatar
Matthieu Constant committed
27
    p.position()[:, 1] += .15/2 
Matthieu Constant's avatar
Matthieu Constant committed
28
    print("Number of grains=%d" % len(p.position()[:, 1]))
Matthieu Constant's avatar
Matthieu Constant committed
29
30
31
32
    p.write_vtk(filename, 0, 0)


outputdir = "output"
Matthieu Constant's avatar
Matthieu Constant committed
33
34
35
36
37
38
39
40
41
if not os.path.isdir(outputdir) :
    os.makedirs(outputdir)

t = 0
ii = 0

p = scontact3.ParticleProblem()

#physical parameters
Matthieu Constant's avatar
Matthieu Constant committed
42
43
44
45
46
g =  -9.81                                      # gravity
rho = 1000                                      # fluid density
rhop = 2640                                     # grains density
nu = 1e-6                                       # kinematic viscosity
tEnd = 100                                      # final time
Matthieu Constant's avatar
Matthieu Constant committed
47
48

#numerical parameters
Matthieu Constant's avatar
Matthieu Constant committed
49
50
51
52
lcmin = 0.001                                   # mesh size
dt = 1e-3                                       # time step
alpha = 2.5e-4                                  # stabilization coefficient
epsilon = alpha*lcmin**2 /nu                    # stabilization parametre
Matthieu Constant's avatar
Matthieu Constant committed
53
54
55
print('epsilon',epsilon)
shutil.copy("mesh.msh", outputdir +"/mesh.msh")

Matthieu Constant's avatar
Matthieu Constant committed
56
#Object particles creation
Matthieu Constant's avatar
Matthieu Constant committed
57
genInitialPosition(outputdir, rhop)
Matthieu Constant's avatar
Matthieu Constant committed
58
p = scontact3.ParticleProblem()
Matthieu Constant's avatar
Matthieu Constant committed
59
p.read_vtk(outputdir,0)
Matthieu Constant's avatar
Matthieu Constant committed
60

Matthieu Constant's avatar
Matthieu Constant committed
61
outf = 100                                      # number of iterations between output files
Matthieu Constant's avatar
Matthieu Constant committed
62

Matthieu Constant's avatar
Matthieu Constant committed
63
64
65
#Object fluid creation + Boundary condition of the fluid (field 0 is horizontal-width velocity; field 1 is vertical velocity; field 2 is the horizontal-depth velocity; field 3 is pressure)
#Format: strong_boundaries = [(Boundary tag, Fluid field, Value)
strong_boundaries = [("Top",0,.0),("Top",1,0.),("Top",2,0.),("Box",1,0.),("Box",0,0.),("Box",2,0.),("Top",3,0.)]
Matthieu Constant's avatar
Matthieu Constant committed
66
fluid = fluid.fluid_problem("mesh.msh",g,nu*rho,rho,epsilon,strong_boundaries)
Matthieu Constant's avatar
Matthieu Constant committed
67
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
Matthieu Constant's avatar
Matthieu Constant committed
68
69

tic = time.clock()
Matthieu Constant's avatar
Matthieu Constant committed
70
fluid.export_vtk(outputdir,0,0)
Matthieu Constant's avatar
Matthieu Constant committed
71
#Computation loop
Matthieu Constant's avatar
Matthieu Constant committed
72
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
Matthieu Constant's avatar
Matthieu Constant committed
73
74
while t < tEnd :
    #Fluid solver
Matthieu Constant's avatar
Matthieu Constant committed
75
76
    fluid.implicit_euler(dt)
    forces = fluid.compute_node_force(dt)
Matthieu Constant's avatar
Matthieu Constant committed
77
    #Computation of the new velocities
Matthieu Constant's avatar
Matthieu Constant committed
78
79
    vn = p.velocity() + forces * dt / p.mass()
    vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
Matthieu Constant's avatar
Matthieu Constant committed
80
    #number of sub time step
Matthieu Constant's avatar
Matthieu Constant committed
81
82
    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()))
Matthieu Constant's avatar
Matthieu Constant committed
83
    #Contact solver
Matthieu Constant's avatar
Matthieu Constant committed
84
    for i in range(nsub) :
Matthieu Constant's avatar
Matthieu Constant committed
85
86
87
        tol = 1e-6
        p.iterate(dt/nsub, forces, tol)
    #Localisation of the grains in the fluid
Matthieu Constant's avatar
Matthieu Constant committed
88
89
    fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
    t += dt
Matthieu Constant's avatar
Matthieu Constant committed
90
    #Output files writting
Matthieu Constant's avatar
Matthieu Constant committed
91
92
93
94
95
96
    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.clock() - tic))