Commit c982e3a9 authored by matthieu's avatar matthieu
Browse files

ajout du maillage et du run correspondant à l'article

parent ea88fc48
L = 0.014;
H = 0.076;
y = 0;
lc = 0.001;
Point(1) = {-L, H, 0, lc};
Point(2) = {-L, -H, 0, lc};
Point(3) = {L, -H, 0, lc};
Point(4) = {L, H, 0, lc};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};
Line Loop(1) = {1:4};
Plane Surface(1) = {1};
Physical Line("Box") = {1,2,3};
Physical Line("Top") = {4};
Physical Surface("Domain") = {1};
Physical Point("PtFix") = {1};
import pyfefp
import scontact2
from pyfefp import scontact2Interface
import numpy as np
import time
import os
from pyfefp.legendre import *
def genInitialPosition(filename, r, lx, ly, rhop) :
p = scontact2.ParticleProblem()
scontact2Interface.MeshLoader(p, "mesh1.msh", ["Top", "Box"])
x = np.arange(lx-r, -lx + r*0.999, -2.1*r)
y = np.arange(ly-r, -ly + r*0.999, -2.1*r)
x, y = np.meshgrid(x, y)
x = x.flat
y = y.flat
for i in range(len(x)) :
p.addParticle((x[i], y[i]), r, r**2 * np.pi * rhop);
p.position()[:, 1] += 0.07542
p.write(filename)
dt = 0.005
tEnd = 110
outputdir = "outputarttestpl1"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
rhop=2500
d = 0.000045
mu = 0.001
rho = 1000
ii = 0
filename = outputdir + "/part-00000"
genInitialPosition(filename, d/2, 0.014, 0.00055, rhop)
## scontact
p = scontact2Interface.ScontactInterface(outputdir+"/part-%05d"%ii)
## lmgc
#from pyfefp import lmgc2Interface
#lmgc2Interface.scontactToLmgc90(filename)
#p = lmgc2Interface.LmgcInterface()
g = [0, -9.81]
print("r = %g, m = %g\n" % (p.r()[0], p.mass()[0]))
print("dtmax = %g Re_p/u = %g" % (p.mass()[0] / (150 * mu * d), d * rho * 0.3/mu))
print("RHOP = %g" % rhop)
outf = 10
fluid = pyfefp.FluidSolver("mesh1.msh", dt, rhop=rhop, rho=rho, g = g, mu=mu, epsilon=1e-04, imex=False, TypeEMu=TriangleP1, TypeEMp=TriangleP1)
fluid.add_boundary_condition("Box", 0, 0.)
fluid.add_boundary_condition("Box", 1, 0.)
fluid.add_boundary_condition("Top", 0, 0.)
fluid.add_boundary_condition("Top", 1, 0.)
fluid.add_boundary_condition("Top", 2, 0.)
fluid.complete()
t = 0
p.write_vtk(outputdir, ii, t)
p.write(outputdir, ii, t)
p.write_boundary_vtk(outputdir, ii, t)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.write_solution(outputdir, ii, t, "vtk")
tic = time.clock()
while t < tEnd :
forces = fluid.solve()
vn = p.velocity() + forces * dt / p.mass()
vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
nsub = max(1, int(np.ceil((vmax * dt * 2)/min(p.r()))))
print("NSUB", nsub,"v", vn[:, 1],"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
#dt=min(0.05,d/vmax)
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, "vtk")
ii += 1
print("%i : %.2g/%.2g (cpu %.2g)" % (ii, t, tEnd, time.clock() - tic))
Markdown is supported
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