depot-fluid.py 2.67 KB
Newer Older
dubois's avatar
maj  
dubois committed
1
2
3
4
5
6
7
8
9
10
11
#!/usr/bin/env python
import pyfefp
from pyfefp.legendre import *
from pyfefp import  scontact2Interface
import scontact2
import numpy as np
from pyfefp import mesh
import os
import time
import random

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
12
def genInitialPosition(filename, r, ox, oy, lx, ly, rhop) :
dubois's avatar
maj  
dubois committed
13
14
    p = scontact2.ParticleProblem()
    scontact2Interface.MeshLoader(p, "mesh.msh", ["Top", "Box"])
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
15
16
    x = np.arange(ox+r, ox+lx, 2*r)
    y = np.arange(oy+r, oy+ly, 2*r)
dubois's avatar
maj  
dubois committed
17
18
19
20
21
22
23
    for i in range(x.shape[0]) :
       for j in range(y.shape[0]) :
          rr=r*(0.8+random.random()*0.2)
          p.addParticle((x[i], y[j]), rr, rr**2 * np.pi * rhop);
    p.write(filename)


Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
24
outputdir = "output-initpressure"
dubois's avatar
maj  
dubois committed
25
26
27
28
29
30
if not os.path.isdir(outputdir) :
    os.makedirs(outputdir)

t = 0
ii = 0

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
31
r=1e-4*3
dubois's avatar
maj  
dubois committed
32
33
34
35
36
37
p = scontact2.ParticleProblem()

#physical parameters
g =  [0, -9.81]
rho = 1200
rhop = 1500 #2400
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
38
nu = 1e-6
dubois's avatar
maj  
dubois committed
39
40
41
tEnd = 20

#numerical parameters
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
42
#h = 2e-4 # approx r*100 but should match the mesh size
dubois's avatar
maj  
dubois committed
43

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
44
c =  20 #h/dt* 7./50.
45
dt = 0.25e-4 #0.5*r/V
dubois's avatar
maj  
dubois committed
46
47
48
49
epsilon = 1e-5 #2e-2*c*h # ?? not sure ??

print("c = %g, eps = %g" % (c, epsilon))

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
50
51
52
53
y0part = 0.051
lypart = 0.025
lxpart = 0.015
genInitialPosition(outputdir + "/part-00000", r, 0., y0part, lxpart, lypart, rhop)
dubois's avatar
maj  
dubois committed
54
55
p = scontact2Interface.ScontactInterface(outputdir+"/part-00000")

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
56
57
58
deltamass = p.volume().sum()*(rhop-rho)
deltapinit = deltamass*-g[1]/lxpart

dubois's avatar
maj  
dubois committed
59
60
print("r = %g, m = %g\n" %  (p.r()[0], p.mass()[0]))
print("RHOP = %g" % rhop)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
61
outf = 100 #000
dubois's avatar
maj  
dubois committed
62

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
63
fluid = pyfefp.FluidSolver("mesh.msh", dt, rhop=rhop, rho=rho, g = g, mu=nu*rho, c=c*c, epsilon = epsilon,timeScheme="Explicit", dragFactor=0.3, TypeEMu=TriangleP1, TypeEMp=TriangleP1)
dubois's avatar
maj  
dubois committed
64
65
66
67
68
69
70
71

# impermeability
fluid.add_boundary_condition("Top", 1, 0.)
fluid.add_boundary_condition("Top", 0, 0.)
fluid.add_boundary_condition("Box", 1, 0.)
fluid.add_boundary_condition("Box", 0, 0.)

fluid.complete()
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
72
73
74
75
76
y = fluid.mesh_position()[...,1].transpose()
pinit = (y0part+lypart-y)*deltapinit/lypart
pinit[pinit<0] = 0
pinit[pinit>deltapinit] = deltapinit
fluid._dof.setField(2, pinit, fluid.solution())
dubois's avatar
maj  
dubois committed
77
78
79
80
81
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())

tic = time.clock()
while t < tEnd :
    forces = fluid.solve() + g * p.mass()
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
82
    NS = 1
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
83
84
    if ii%NS == 0 :
        p.iterate(dt*NS, forces)
dubois's avatar
maj  
dubois committed
85
86
87
88
    fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
    t += dt
    if ii %outf == 0 :
        ioutput = int(ii/outf) + 1
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
89
        #p.write_vtk(outputdir, ioutput, t)
dubois's avatar
maj  
dubois committed
90
        p.write(outputdir, ioutput, t)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
91
92
93
        #p.write_boundary_vtk(outputdir, ioutput, t)
        fluid.write_solution(outputdir, ioutput, t, "msh")
        #fluid.write_solution(outputdir, ioutput, t, "vtk")
dubois's avatar
maj  
dubois committed
94
95
96
    ii += 1
    print("%i : %.2g/%.2g (cpu %.2g)" % (ii, t, tEnd, time.clock() - tic))