unsortedBidisperse.py 3.84 KB
Newer Older
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
def genInitialPosition(filename, r1, r2, rout, rhop1, rhop2) :
12
13
    p = scontact2.ParticleProblem()
    p.load_msh_boundaries("mesh.msh", ["Top", "Bottom","Left","Right"])
Matthieu Constant's avatar
Matthieu Constant committed
14
15
16
    e = 1.4e-3/sqrt(2)
    x = np.arange(rout, -rout, -1.4e-3)
    y = np.arange(rout, -rout, -1.4e-3)
17
18
    x, y = np.meshgrid(x, y)
    R2 = x**2 + y**2
Matthieu Constant's avatar
Matthieu Constant committed
19
    keep = R2 < (rout - r2/2)**2
20
21
22
23
24
25
    x = x[keep]
    y = y[keep]
    N1 = 0
    N2 = 0
    for i in range(x.shape[0]) :
        a = random.choice([1,2])
Matthieu Constant's avatar
Matthieu Constant committed
26
        if a == 1 and N1 <= 17:
27
            N1 += 1
Matthieu Constant's avatar
Matthieu Constant committed
28
29
            p.add_particle((x[i]+random.uniform(-1.4e-3/3,1.4e-3/3), y[i]+random.uniform(-1.4e-3/3,1.4e-3/3)), r1, r1**2 * np.pi * rhop1);
        elif a == 2 and N2 <= 17:
30
            N2 += 1
Matthieu Constant's avatar
Matthieu Constant committed
31
            p.add_particle((x[i]+random.uniform(-1.4e-3/3,1.4e-3/3), y[i]+random.uniform(-1.4e-3/3,1.4e-3/3)), r2, r2**2 * np.pi * rhop2);
32
        else:
Matthieu Constant's avatar
Matthieu Constant committed
33
34
35
36
            if N1 > 17 and a == 1:
                p.add_particle((x[i]+random.uniform(-1.4e-3/3,1.4e-3/3), y[i]+random.uniform(-1.4e-3/3,1.4e-3/3)), r2, r2**2 * np.pi * rhop2);
            if N2 > 17 and a == 2:
                p.add_particle((x[i]+random.uniform(-1.4e-3/3,1.4e-3/3), y[i]+random.uniform(-1.4e-3/3,1.4e-3/3)), r1, r1**2 * np.pi * rhop1);
37
38
39
40
41
42
43
44

        
    print('N1',N1,'N2',N2,'N',N1+N2)
    p.position()[:, 1] += 0.45
    p.write(filename)



Matthieu Constant's avatar
Matthieu Constant committed
45
outputdir = "outputGoute"
46
47
48
49
50
51
if not os.path.isdir(outputdir) :
    os.makedirs(outputdir)

t = 0
ii = 0

Matthieu Constant's avatar
Matthieu Constant committed
52
53
54
r1=0.36e-3/2/2
r2 = 0.78e-3/2/2
R = 4.1e-3/2
55
56
57
58
59
60
p = scontact2.ParticleProblem()
#R = np.random.uniform(45e-06, 90e-06, len(x))

#physical parameters
g =  -9.81
rho = 950
Matthieu Constant's avatar
Matthieu Constant committed
61
62
rhop1 = 2440
rhop2 = 2550
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
nu = 174e-6#/(10**0.5)
V = 0.5 # todo : estimate V base on limit velocity
print('V',V)
tEnd = 10000

#numerical parameters
lcmin = 0.005 # approx r*100 but should match the mesh size
dt = 1e-2
alpha = 2.5e-4
epsilon = alpha*lcmin**2 /nu#2e-2*c*h*2#2e-2*c*h # ?? not sure ??1e-5
autoEpsilon = False
print('epsilon',epsilon)

shutil.copy("mesh.msh", outputdir +"/mesh.msh")
#scontact2Interface.MeshLoader(p, "funnel.msh", ("Funnel", "Top", "Bottom", "Lateral"))
p.write(outputdir+"/part-00000")

Matthieu Constant's avatar
Matthieu Constant committed
80
genInitialPosition(outputdir + "/part-00000", r1, r2, R, rhop1, rhop2)
81
82
83
84

p = scontact2.ParticleProblem(outputdir+"/part-00000")

print("r = %g, m = %g\n" %  (p.r()[0], p.mass()[0]))
Matthieu Constant's avatar
Matthieu Constant committed
85
print("RHOP = %g" % rhop1)
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
outf = 1
outf1 = 100000

strong_boundaries = [("Top",2,0.),("Top",1,0.),("Top",0,0.),("Bottom",1,0.),("Bottom",0,0.),("Left",0,0.),("Left",1,0.),("Right",0,0.),("Right",1,0.)]
fluid = fluid.fluid_problem("mesh.msh",g,nu*rho,rho,epsilon,autoEpsilon,strong_boundaries)

fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.export(outputdir,0,0)

ii = 0
    

tic = time.clock()
forces = g*p.mass()
while t < tEnd : 
   # if (ii%50==0 and ii != 0):
    #   fluid.adapt_mesh(0.01,50,4e-5,autoEpsilon)
Matthieu Constant's avatar
Matthieu Constant committed
103
    forces = fluid.compute_node_force(dt,1e-3)
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
    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(10*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)  
    fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
    fluid.implicit_euler(dt)
    t += dt
    if ii %outf == 0 :
        ioutput = int(ii/outf) + 1
        p.write_vtk(outputdir, ioutput, t)
        p.write(outputdir+"/part-%05d" % ioutput)
        p.write_boundary_vtk(outputdir, ioutput, t)
        fluid.export(outputdir, t, ioutput)
    ii += 1
    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.clock() - tic))