drop.py 2.8 KB
Newer Older
Matthieu Constant's avatar
Matthieu Constant committed
1
#!/usr/bin/env python
2
from marblesbag import fluid as fluid
3
4
from marblesbag import scontact2

Matthieu Constant's avatar
Matthieu Constant committed
5
6
7
8
import numpy as np
import os
import time
import shutil
9
import random
Matthieu Constant's avatar
Matthieu Constant committed
10
11
12

def genInitialPosition(filename, r, rout, rhop) :
    p = scontact2.ParticleProblem()
13
14
    p.load_msh_boundaries("mesh.msh", ["Top", "Box"])
    x = np.arange(rout, -rout, -2*r*2**0.5)
Matthieu Constant's avatar
Matthieu Constant committed
15
16
17
18
19
20
    x, y = np.meshgrid(x, x)
    R2 = x**2 + y**2
    keep = R2 < (rout - r)**2
    x = x[keep]
    y = y[keep]
    for i in range(x.shape[0]) :
21
22
23
        p.add_particle((x[i], y[i]), r, r**2 * np.pi * rhop);
    
    p.position()[:, 1] += 0.45
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
24
    p.write_vtk(filename,0,0)
Matthieu Constant's avatar
Matthieu Constant committed
25
26


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

t = 0
ii = 0

34
r=25e-6/2
Matthieu Constant's avatar
Matthieu Constant committed
35
36
37
38
p = scontact2.ParticleProblem()
#R = np.random.uniform(45e-06, 90e-06, len(x))

#physical parameters
39
g =  -9.81
Matthieu Constant's avatar
Matthieu Constant committed
40
41
rho = 1200
rhop = 2400
42
nu = 0.8e-4
Matthieu Constant's avatar
Matthieu Constant committed
43
44
V = 0.5 # todo : estimate V base on limit velocity
print('V',V)
45
tEnd = 100
Matthieu Constant's avatar
Matthieu Constant committed
46
47

#numerical parameters
48
lcmin = 0.0001 # approx r*100 but should match the mesh size
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
49
dt = 5e-4*10
50
alpha = 2.5e-2
Matthieu Constant's avatar
drop2d    
Matthieu Constant committed
51
epsilon = alpha*lcmin**2 /nu
52
print('epsilon',epsilon)
Matthieu Constant's avatar
Matthieu Constant committed
53
54
55
56

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

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
57
genInitialPosition(outputdir, r, 1e-3, rhop)
Matthieu Constant's avatar
Matthieu Constant committed
58

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
59
60
p = scontact2.ParticleProblem()
p.read_vtk(outputdir,0)
Matthieu Constant's avatar
Matthieu Constant committed
61
62
63

print("r = %g, m = %g\n" %  (p.r()[0], p.mass()[0]))
print("RHOP = %g" % rhop)
Matthieu Constant's avatar
drop2d    
Matthieu Constant committed
64
outf = 10
Matthieu Constant's avatar
Matthieu Constant committed
65
66
67
outf1 = 100000

strong_boundaries = [("Top",2,0.),("Top",1,0.),("Top",0,0.),("Box",1,0.),("Box",0,0.)]
Matthieu Constant's avatar
drop2d    
Matthieu Constant committed
68
fluid = fluid.fluid_problem("mesh.msh",g,nu*rho,rho,epsilon,strong_boundaries)
Matthieu Constant's avatar
Matthieu Constant committed
69
70

fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
71
fluid.export_vtk(outputdir,0,0)
Matthieu Constant's avatar
Matthieu Constant committed
72

73
ii = 0
Matthieu Constant's avatar
Matthieu Constant committed
74
75
tic = time.clock()
forces = g*p.mass()
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
76
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
Matthieu Constant's avatar
Matthieu Constant committed
77
while t < tEnd : 
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
78
    if (ii%5==0 and ii != 0):
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
79
       fluid.adapt_mesh(1e-3,1e-4,10000)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
80
81
82
83
       fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
    else :
       fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
    fluid.implicit_euler(dt)
Matthieu Constant's avatar
Matthieu Constant committed
84
85
86
87
88
89
90
91
92
93
94
    forces = fluid.compute_node_force(dt)
    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)  
    t += dt
    if ii %outf == 0 :
        ioutput = int(ii/outf) + 1
        p.write_vtk(outputdir, ioutput, t)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
95
        #p.write(outputdir+"/part-%05d" % ioutput)
96
        fluid.export_vtk(outputdir, t, ioutput)
Matthieu Constant's avatar
Matthieu Constant committed
97
98
    ii += 1
    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.clock() - tic))