Commit 64b384e5 authored by dubois's avatar dubois
Browse files

maj

parent adf2a82d
#!/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()
d=1e-4;
d=1e-3;
L = 15*d;
H = 50*d;
lc = d;
lc = 2*d;
R=3*d;
Point(1) = {0., 0., 0, lc};
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment