drop.py 3.76 KB
Newer Older
Matthieu Constant's avatar
Matthieu Constant committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
#!/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

def genInitialPosition(filename, r1, r2, rout, rhop1, rhop2) :
    p = scontact2.ParticleProblem()
    p.load_msh_boundaries("mesh.msh", ["Top", "Bottom","Left","Right"])
    e = .7e-4
    x = np.arange(rout, -rout, -e)
    y = np.arange(rout, -rout, -e)
    x, y = np.meshgrid(x, y)
    R2 = x**2 + y**2
    keep = R2 < (rout - r2/2)**2
    x = x[keep]
    y = y[keep]
    N1 = 0
    N2 = 0
    for i in range(x.shape[0]) :
        a = random.choice([1,2])
        if a == 1 and N1 <= 128:
            N1 += 1
            p.add_particle((x[i]+random.uniform(-e/3,e/3), y[i]+random.uniform(-e/3,e/3)), r1, r1**2 * np.pi * rhop1);
        elif a == 2 and N2 <= 128:
            N2 += 1
            p.add_particle((x[i]+random.uniform(-e/3,e/3), y[i]+random.uniform(-e/3,e/3)), r2, r2**2 * np.pi * rhop2);
        else:
            if N1 > 128 and a == 1:
                p.add_particle((x[i]+random.uniform(-e/3,e/3), y[i]+random.uniform(-e/3,e/3)), r2, r2**2 * np.pi * rhop2);
            if N2 > 128 and a == 2:
                p.add_particle((x[i]+random.uniform(-e/3,e/3), y[i]+random.uniform(-e/3,e/3)), r1, r1**2 * np.pi * rhop1);

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


Matthieu Constant's avatar
Matthieu Constant committed
44
outputdir = "outputGoutteStokesGoodNew"
Matthieu Constant's avatar
Matthieu Constant committed
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
if not os.path.isdir(outputdir) :
    os.makedirs(outputdir)

t = 0
ii = 0

r1=25e-6/2
r2 = 25e-6/4
R = 1e-3
p = scontact2.ParticleProblem()
#R = np.random.uniform(45e-06, 90e-06, len(x))

#physical parameters
g =  -9.81
rho = 1200
rhop1 = 2400
rhop2 = 2400
Matthieu Constant's avatar
Matthieu Constant committed
62
nu = 1e-4#/(100**0.5)
Matthieu Constant's avatar
Matthieu Constant committed
63
64
V = 0.5 # todo : estimate V base on limit velocity
print('V',V)
Matthieu Constant's avatar
Matthieu Constant committed
65
tEnd = 10
Matthieu Constant's avatar
Matthieu Constant committed
66
67
68

#numerical parameters
lcmin = 0.00005 # approx r*100 but should match the mesh size
Matthieu Constant's avatar
Matthieu Constant committed
69
dt = 5e-4
Matthieu Constant's avatar
Matthieu Constant committed
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
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")

genInitialPosition(outputdir + "/part-00000", r1, r2, R, rhop1, rhop2)

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

print("r = %g, m = %g\n" %  (p.r()[0], p.mass()[0]))
print("RHOP = %g" % rhop1)
Matthieu Constant's avatar
Matthieu Constant committed
85
outf = 50
Matthieu Constant's avatar
Matthieu Constant committed
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
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 : 
Matthieu Constant's avatar
Matthieu Constant committed
101
102
103
    if (ii%100==0 and ii != 0) or ii==10:
       fluid.adapt_mesh(0.01,100,5e-5,autoEpsilon)
    forces = fluid.compute_node_force(dt,10*r1)
Matthieu Constant's avatar
Matthieu Constant committed
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))