fd_depot-fluid.py 2.69 KB
Newer Older
Frédéric Dubois's avatar
en vrac  
Frédéric Dubois committed
1
2
#!/usr/bin/env python

Frédéric Dubois's avatar
Frédéric Dubois committed
3
use_lmgc90 = False #True
Frédéric Dubois's avatar
en vrac  
Frédéric Dubois committed
4

Matthieu Constant's avatar
Matthieu Constant committed
5
6
from migflow import fluid
from migflow import scontact
Frédéric Dubois's avatar
Frédéric Dubois committed
7
8
if use_lmgc90 :
    from migflow import lmgc90Interface
Frédéric Dubois's avatar
en vrac  
Frédéric Dubois committed
9
10
11
12
13
14
15
16
17
18
19
20

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

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

def genInitialPosition(filename, r, ox, oy, lx, ly, rhop) :
Matthieu Constant's avatar
Matthieu Constant committed
21
    p = scontact.ParticleProblem(2)
Frédéric Dubois's avatar
en vrac  
Frédéric Dubois committed
22
23
24
25
26
27
28
    p.load_msh_boundaries("mesh.msh", ["Left", "Box"])
    x = np.arange(ox+r, ox+lx, 2*r)
    y = np.arange(oy+r, oy+ly, 2*r)
    for i in range(x.shape[0]) :
       for j in range(y.shape[0]) :
          rr=r*(0.95+random.random()*0.05)
          p.add_particle((x[i], y[j]), rr, rr**2 * np.pi * rhop);
Matthieu Constant's avatar
Matthieu Constant committed
29
    p.write_vtk(filename,0,0)
Frédéric Dubois's avatar
en vrac  
Frédéric Dubois committed
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51


t = 0
ii = 0

#physical parameters
g =  0. #-9.81
rho = 1200
rhop = 1800 #2400
nu = 1e-4
tEnd = 20.

#numerical parameters
#h = 2e-4 # approx r*100 but should match the mesh size

dt = 1e-3      #0.5*r/V

r=1e-4
x0part = 0
y0part = 0
lypart = 150*r
lxpart = 50*r
Matthieu Constant's avatar
Matthieu Constant committed
52
genInitialPosition(outputdir, r, x0part, y0part, lxpart, lypart, rhop)
Frédéric Dubois's avatar
en vrac  
Frédéric Dubois committed
53

Frédéric Dubois's avatar
Frédéric Dubois committed
54
55
56
57
58
if use_lmgc90 :
    # friction coefficient    
    friction = 0.1                                    
    lmgc90Interface.scontactTolmgc90(outputdir,2,0,friction)
    p = lmgc90Interface.ParticleProblem(2)
Frédéric Dubois's avatar
en vrac  
Frédéric Dubois committed
59
else :
Matthieu Constant's avatar
Matthieu Constant committed
60
61
    p = scontact.ParticleProblem(2)
    p.read_vtk(outputdir,0)
Frédéric Dubois's avatar
en vrac  
Frédéric Dubois committed
62
63
64
65
66
67
68
69

deltamass = p.volume().sum()*(rhop-rho)
deltapinit = deltamass*-g/lxpart

print("r = %g, m = %g\n" %  (p.r()[0], p.mass()[0]))
print("RHOP = %g" % rhop)
outf = 20 #000

Matthieu Constant's avatar
Matthieu Constant committed
70
71
72
73
74
75
76
fluid = fluid.FluidProblem(2,g,nu*rho,rho)
fluid.load_msh("mesh.msh")
fluid.set_strong_boundary("Box",0,0)
fluid.set_strong_boundary("Box",1,0)
fluid.set_strong_boundary("Left",0,0.001)
fluid.set_strong_boundary("Right",2,0)
fluid.set_weak_boundary("Left","Inflow")
Frédéric Dubois's avatar
Frédéric Dubois committed
77
fluid.set_weak_boundary("Right","Outflow")
Frédéric Dubois's avatar
en vrac  
Frédéric Dubois committed
78
79
80
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.export_vtk(outputdir,0,0)

Matthieu Constant's avatar
Matthieu Constant committed
81
tic = time.time()
Frédéric Dubois's avatar
en vrac  
Frédéric Dubois committed
82
83
84
while t < tEnd :
    #if ii%100==0 and ii != 0:
     #  fluid.adapt_mesh(0.1,10,8e-4)
Matthieu Constant's avatar
Matthieu Constant committed
85
    fluid.implicit_euler(dt)
Frédéric Dubois's avatar
en vrac  
Frédéric Dubois committed
86
87
88
89
90
91
92
93
94
95
96
    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)  
    fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
    t += dt
    if ii %outf == 0 :
        ioutput = int(ii/outf) + 1
Frédéric Dubois's avatar
Frédéric Dubois committed
97
        if not use_lmgc90 : p.write_vtk(outputdir, ioutput, t)
Frédéric Dubois's avatar
en vrac  
Frédéric Dubois committed
98
99
        fluid.export_vtk(outputdir, t, ioutput)
    ii += 1
Matthieu Constant's avatar
Matthieu Constant committed
100
    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
Frédéric Dubois's avatar
en vrac  
Frédéric Dubois committed
101