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

import numpy as np
import os
import time
import shutil
import random

Matthieu Constant's avatar
Matthieu Constant committed
11
12
13
14
#grains creation and placement + physical boundaries definition
'''r is the radius of the grains
ly is the height of the box
rhop is the particles density'''
Matthieu Constant's avatar
Matthieu Constant committed
15
16
def genInitialPosition(filename, r, ly, rhop) :
    p = scontact2.ParticleProblem()
Matthieu Constant's avatar
Matthieu Constant committed
17
    #Loading of the mesh.msh file specifying physical boundaries name
Matthieu Constant's avatar
Matthieu Constant committed
18
    p.load_msh_boundaries("mesh.msh", ["Top", "Lateral","Bottom","BottomOut","TopOut"])
Matthieu Constant's avatar
Matthieu Constant committed
19
    #Definition of the points where the grains are located
Matthieu Constant's avatar
Matthieu Constant committed
20
21
22
23
24
25
    x = np.arange(-0.025+r, 0.025-r, r)
    y = np.arange(0.06-r, .04, -r)
    x, y = np.meshgrid(x, y)
    x = x.flat
    y = y.flat
    for i in range(len(x)) :
Matthieu Constant's avatar
Matthieu Constant committed
26
        rhop = random.choice([1100,1200,1300])
Matthieu Constant's avatar
Matthieu Constant committed
27
        #Addition of an particle object at each point
Matthieu Constant's avatar
Matthieu Constant committed
28
        p.add_particle((x[i], y[i]), r/2, r**2/4 * np.pi * rhop);
Matthieu Constant's avatar
Matthieu Constant committed
29
    p.write_vtk(filename,0,0)
Matthieu Constant's avatar
Matthieu Constant committed
30
31


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

t = 0
ii = 0


#physical parameters
Matthieu Constant's avatar
Matthieu Constant committed
41
42
43
44
45
46
47
r = 2e-4                            # radius
ly = 5e-2                           # box height
g =  -9.81                          # gravity
rho = 1000                          # fluid density
rhop = 1500                         # grains density
nu = 1e-6                           # kinematic viscosity
tEnd = 100                          # final time
Matthieu Constant's avatar
Matthieu Constant committed
48
49

#numerical parameters
Matthieu Constant's avatar
Matthieu Constant committed
50
51
52
53
lcmin = 0.001                       # mesh size
dt = 1e-2                           # time step
alpha = 2.5e-5                      # stabilization coefficient
epsilon = alpha*lcmin**2 /nu        # stabilization parameter
Matthieu Constant's avatar
Matthieu Constant committed
54
55
56
57
print('epsilon',epsilon)

shutil.copy("mesh.msh", outputdir +"/mesh.msh")

Matthieu Constant's avatar
Matthieu Constant committed
58
#Object particles creation
Matthieu Constant's avatar
Matthieu Constant committed
59
60
61
genInitialPosition(outputdir, r, ly, rhop)
p = scontact2.ParticleProblem()
p.read_vtk(outputdir,0)
Matthieu Constant's avatar
Matthieu Constant committed
62
63
64

print("r = %g, m = %g\n" %  (p.r()[0], p.mass()[0]))
print("RHOP = %g" % rhop)
Matthieu Constant's avatar
Matthieu Constant committed
65
outf = 20                           # number of iterations between output files
Matthieu Constant's avatar
Matthieu Constant committed
66
67
68

ii = 0

Matthieu Constant's avatar
Matthieu Constant committed
69
70
#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)
Matthieu Constant's avatar
Matthieu Constant committed
71
72
73
74
75
76
77
strong_boundaries = [("Top",2,0.),("TopOut",1,0),("Top",1,0),("BottomOut",1,0),("Bottom",1,0),("Lateral",0,0.),("Lateral",1,0.)]
fluid = fluid.fluid_problem("mesh.msh",g,nu*rho,rho,epsilon,strong_boundaries)
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
78
#Computation loop
Matthieu Constant's avatar
Matthieu Constant committed
79
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
Matthieu Constant's avatar
Matthieu Constant committed
80
81
while t < tEnd :
    #Fluid solver
Matthieu Constant's avatar
Matthieu Constant committed
82
    fluid.implicit_euler(dt)
Matthieu Constant's avatar
Matthieu Constant committed
83
    forces = fluid.compute_node_force(dt)
Matthieu Constant's avatar
Matthieu Constant committed
84
    #Computation of the new velocities
Matthieu Constant's avatar
Matthieu Constant committed
85
86
    vn = p.velocity() + forces * dt / p.mass()
    vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
Matthieu Constant's avatar
Matthieu Constant committed
87
    #number of sub time step
Matthieu Constant's avatar
Matthieu Constant committed
88
89
    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
90
    #Contact solver
Matthieu Constant's avatar
Matthieu Constant committed
91
    for i in range(nsub) :
Matthieu Constant's avatar
Matthieu Constant committed
92
93
94
        tol = 1e-6
        p.iterate(dt/nsub, forces, tol)
    #Localisation of the grains in the fluid
Matthieu Constant's avatar
Matthieu Constant committed
95
96
    fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
    t += dt
Matthieu Constant's avatar
Matthieu Constant committed
97
    #Output files writting
Matthieu Constant's avatar
Matthieu Constant committed
98
99
100
101
102
103
    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))