depot.py 1.94 KB
Newer Older
Matthieu Constant's avatar
depot    
Matthieu Constant committed
1
#!/usr/bin/env python
Matthieu Constant's avatar
Matthieu Constant committed
2
3
from migflow import fluid
from migflow import scontact
Matthieu Constant's avatar
depot    
Matthieu Constant committed
4
5
6
7
8
9
10
11

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

def genInitialPosition(filename, r, N, rhop) :
Matthieu Constant's avatar
Matthieu Constant committed
12
    p = scontact.ParticleProblem(2)
Matthieu Constant's avatar
depot    
Matthieu Constant committed
13
14
    p.load_msh_boundaries("mesh.msh", ["Top", "Bottom","Lateral","Injection"])
    dr = r/10
Matthieu Constant's avatar
fix dep    
Matthieu Constant committed
15
16
    x = np.arange(r+dr, .222-r-dr, 2*(r+dr))
    y = np.arange(r+dr, .4-r-dr, 2*(r+dr))
Matthieu Constant's avatar
depot    
Matthieu Constant committed
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
    i = 0
    for yl in y :
        for xl in x:
            if i > N:
                break
            r1 = random.choice([0.98*r,r,1.02*r])
            d=np.random.random()
            p.add_particle((xl+(d*dr), yl+(d*dr)), r1, r1**2 * np.pi * rhop);
            i += 1
    print(len(p.mass()))
    p.write_vtk(filename,0,0)


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

t = 0
ii = 0

r = 5e-4
ly = 5e-2
Matthieu Constant's avatar
Matthieu Constant committed
39
p = scontact.ParticleProblem(2)
Matthieu Constant's avatar
depot    
Matthieu Constant committed
40

Matthieu Constant's avatar
Matthieu Constant committed
41
# physical parameters
Matthieu Constant's avatar
depot    
Matthieu Constant committed
42
43
44
45
46
47
g =  -9.81
rho = 1000
rhop = 2500
nu = 1e-6
tEnd = 2.5

Matthieu Constant's avatar
Matthieu Constant committed
48
# numerical parameters
Matthieu Constant's avatar
depot    
Matthieu Constant committed
49
50
51
52
dt = 1e-4
N = 40000
genInitialPosition(outputdir, r, N, rhop)

Matthieu Constant's avatar
Matthieu Constant committed
53
p = scontact.ParticleProblem(2)
Matthieu Constant's avatar
depot    
Matthieu Constant committed
54
55
56
57
58
59
60
61
62
p.read_vtk(outputdir,0)

print("r = %g, m = %g\n" %  (p.r()[0], p.mass()[0]))
print("RHOP = %g" % rhop)
outf = 200
outf1 = 100000

ii = 0

63
forces = np.zeros((p.n_particles(),p.dim()))
Matthieu Constant's avatar
depot    
Matthieu Constant committed
64
print(forces[:,1].shape)
65
tic = time.process_time()
Matthieu Constant's avatar
depot    
Matthieu Constant committed
66
67
68
69
70
71
72
73
74
75
76
77
78
while t < tEnd : 
    forces[:,1] = g*p.mass()[0]-p.velocity()[:,1]
    vn = p.velocity() + forces * dt / p.mass() 
    vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
    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()))
    for i in range(nsub) :
        p.iterate(dt/nsub, forces) 
    t += dt
    if ii %outf == 0 :
        ioutput = int(ii/outf) + 1
        p.write_vtk(outputdir, ioutput, t)
    ii += 1
79
    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.process_time() - tic))
Matthieu Constant's avatar
depot    
Matthieu Constant committed
80
81