depot-fluid.py 3.74 KB
Newer Older
dubois's avatar
maj  
dubois committed
1
#!/usr/bin/env python
2
3
from marblesbag import fluid as fluid
from marblesbag import scontact2
Matthieu Constant's avatar
Matthieu Constant committed
4
from marblesbag import lmgc2Interface
5

dubois's avatar
maj  
dubois committed
6
7
8
import numpy as np
import os
import time
9
import shutil
dubois's avatar
maj  
dubois committed
10
import random
Matthieu Constant's avatar
Matthieu Constant committed
11
#grains creation and placement + physical boundaries definition
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
12
def genInitialPosition(filename, r, ox, oy, lx, ly, rhop) :
dubois's avatar
maj  
dubois committed
13
    p = scontact2.ParticleProblem()
Matthieu Constant's avatar
Matthieu Constant committed
14
    #Loading of the mesh.msh file specifying physical boundaries name
15
    p.load_msh_boundaries("mesh.msh", ["Top", "Box"])
Matthieu Constant's avatar
Matthieu Constant committed
16
    #Definition of the points where the grains are located
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
17
18
    x = np.arange(ox+r, ox+lx, 2*r)
    y = np.arange(oy+r, oy+ly, 2*r)
dubois's avatar
maj  
dubois committed
19
20
    for i in range(x.shape[0]) :
       for j in range(y.shape[0]) :
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
21
          rr=r*(0.95+random.random()*0.05)
Matthieu Constant's avatar
Matthieu Constant committed
22
          #Addition of an particle object at each point
23
          p.add_particle((x[i], y[j]), rr, rr**2 * np.pi * rhop);
Matthieu Constant's avatar
Matthieu Constant committed
24
    p.write_vtk(filename,0,0)
dubois's avatar
maj  
dubois committed
25
26


27
outputdir = "output"
dubois's avatar
maj  
dubois committed
28
29
30
31
32
33
34
if not os.path.isdir(outputdir) :
    os.makedirs(outputdir)

t = 0
ii = 0

#physical parameters
Matthieu Constant's avatar
Matthieu Constant committed
35
36
37
38
39
40
41
42
43
g =  -9.81                                      # gravity
rho = 1200                                      # fluid density
rhop = 1800                                     # grains density
nu = 1e-6                                       # kinematic viscosity
tEnd = 7                                        # final time
r=1e-4                                          # grains radius
y0part = 40e-3                                  # vertical center of the grains
lypart = 20e-3                                  # semi-height of the grains
lxpart = 0.015                                  # semi-width of the grains
dubois's avatar
maj  
dubois committed
44
45

#numerical parameters
Matthieu Constant's avatar
Matthieu Constant committed
46
47
dt = 1e-3                                       # time step
epsilon = 1e-4                                  # stabilization parameter
dubois's avatar
maj  
dubois committed
48

Matthieu Constant's avatar
Matthieu Constant committed
49
print("eps = %g" % ( epsilon))
Matthieu Constant's avatar
Matthieu Constant committed
50
#Enable the use of lmgc90
Matthieu Constant's avatar
Matthieu Constant committed
51
use_lmgc = True
Matthieu Constant's avatar
Matthieu Constant committed
52
53
54


#Object particles creation
Matthieu Constant's avatar
Matthieu Constant committed
55
genInitialPosition(outputdir, r, 0., y0part, lxpart, lypart, rhop)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
56
if use_lmgc :
Matthieu Constant's avatar
Matthieu Constant committed
57
58
    friction=0.1                                # friction coefficient
    lmgc2Interface.scontactToLmgc90(outputdir, 0, friction)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
59
60
    p = lmgc2Interface.LmgcInterface()
else :
Matthieu Constant's avatar
Matthieu Constant committed
61
62
    p = scontact2.ParticleProblem()
    p.read_vtk(outputdir,0)
dubois's avatar
maj  
dubois committed
63
64

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


Matthieu Constant's avatar
Matthieu Constant committed
68
69
#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)
70
strong_boundaries = [("Box",0,0),("Box",1,0),("Top",0,0.),("Top",1,0.),("Top",2,0.)]
Matthieu Constant's avatar
Matthieu Constant committed
71
fluid = fluid.fluid_problem("mesh.msh",g,nu*rho,rho,epsilon,strong_boundaries)
dubois's avatar
maj  
dubois committed
72
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
73
fluid.export(outputdir,0,0)
dubois's avatar
maj  
dubois committed
74
75

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