depot-fluid.py 3.31 KB
Newer Older
dubois's avatar
maj  
dubois 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
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
#!/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

def genInitialPosition(filename, r, ox, oy, nx, ny, rhop) :
    p = scontact2.ParticleProblem()
    scontact2Interface.MeshLoader(p, "mesh.msh", ["Top", "Box"])
    x = np.arange(ox+2*r, ox+(2*nx*r), 2*r)
    y = np.arange(oy+2*r, oy+2*ny*r, 2*r)

    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)


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

t = 0
ii = 0

d=1e-3 # la valeur de ref dans le .geo
r=0.1*d
p = scontact2.ParticleProblem()

#physical parameters
g =  [0, -9.81]
rho = 1200
rhop = 1500 #2400
nu = 1e-4
V = 1e-3  #0.5 # todo : estimate V base on limit velocity
print('V',V)
tEnd = 20

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

dt = 5e-4 #0.5*r/V

c =  1e-2 #h/dt* 7./50.

epsilon = 1e-5 #2e-2*c*h # ?? not sure ??

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

genInitialPosition(outputdir + "/part-00000", r, 0., 51*d, 75, 200, rhop)

p = scontact2Interface.ScontactInterface(outputdir+"/part-00000")

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

fluid = pyfefp.FluidSolver("mesh.msh", dt, rhop=rhop, rho=rho, g = g, mu=nu*rho, c=c*c, epsilon = epsilon,timeScheme="Explicit", TypeEMu=TriangleP1, TypeEMp=TriangleP1)

# 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()

fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())

filenergy =open(outputdir+'/filenergy.txt','w')
filenergy.write('Energie cinetique solide, Energie potentielle solide, Energie cinetique fluide, Energie totale\n')

p.write_vtk(outputdir, 1, 0.)
p.write_boundary_vtk(outputdir, 1, 0.)

tic = time.clock()
while t < tEnd :
    forces = fluid.solve() + g * p.mass()
    vn = p.velocity() + forces * dt / p.mass()

    vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
    Ecs = sum(np.hypot(vn[:, 0], vn[:, 1])**2 *0.5)
    Eps = sum(rhop*p.position()[:,1]*9.81)
    if ii % outf1 == 0:
        filenergy.write('%.12f, %.12f, %.12f, %.12f\n'% (Ecs,Eps,fluid._law.Ecf,Ecs+Eps+fluid._law.Ecf))
    if (vmax > 10*V) :
        print("CRASH !!")
        #exit()
    if ii==10000:
        fluid._law.epsilon=2e-1*c*h
    #print("NSUB", nsub,
    print("VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))        
    nsub = max(1, int(np.ceil((vmax * dt * 4)/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
    print("dt",dt)
    if ii %outf == 0 :
        ioutput = int(ii/outf) + 1
        p.write_vtk(outputdir, ioutput, t)
        p.write(outputdir, ioutput, t)
        p.write_boundary_vtk(outputdir, ioutput, t)
        #fluid.write_solution(outputdir, ioutput, t, "msh")
        fluid.write_solution(outputdir, ioutput, t, "vtk")
    ii += 1
    print("%i : %.2g/%.2g (cpu %.2g)" % (ii, t, tEnd, time.clock() - tic))

filenergy.close()