Commit a1e49986 authored by Matthieu Constant's avatar Matthieu Constant
Browse files

wip debut depot 3d

parent 02db5e4f
Pipeline #4542 passed with stage
in 2 minutes and 3 seconds
......@@ -38,6 +38,7 @@ def genInitialPosition(filename, rM, rm, ly, lxz, rhop) :
for z in np.arange(-lxz+rM, lxz-rM, 2*rM):
r = np.random.uniform(rm,rM)
p.add_particle((x, y, z), r, 4.* r**3 * np.pi * rhop / 3.);
print("nb of grains = %d"%len(p.mass()))
p.write_vtk(filename,0,0)
......@@ -49,7 +50,7 @@ if not os.path.isdir(outputdir) :
g = -9.81
rho = 1000
rhop = 2500
nu = 1e-8
nu = 1e-6
V = 0.5 # todo : estimate V base on limit velocity
print('V',V)
tEnd = 5.
......@@ -66,34 +67,30 @@ lxz = 0.025
genInitialPosition(outputdir, rM, rm, ly, lxz, rhop)
friction = 0.1
lmgc90Interface.scontactTolmgc90(outputdir,3,0,friction,assume_box=True)
p = lmgc90Interface.ParticleProblem(3)
#lmgc90Interface.scontactTolmgc90(outputdir,3,0,friction,assume_box=True)
#p = lmgc90Interface.ParticleProblem(3)
p = scontact.ParticleProblem(3)
p.read_vtk(outputdir,0)
#numerical parameters
dt = 1e-4 # time step
fluid = fluid.FluidProblem(3,g,[nu*rho],[rho])
fluid.load_msh("mesh.msh")
fluid.set_strong_boundary("Top",0,0.)
fluid.set_strong_boundary("Top",1,0.)
fluid.set_strong_boundary("Top",2,0.)
fluid.set_strong_boundary("X",1,0.)
fluid.set_strong_boundary("X",2,0.)
fluid.set_strong_boundary("X",0,0.)
fluid.set_strong_boundary("Z",1,0.)
fluid.set_strong_boundary("Z",2,0.)
fluid.set_strong_boundary("Z",0,0.)
fluid.set_strong_boundary("Bottom",0,0.)
fluid.set_strong_boundary("Bottom",1,0.)
fluid.set_strong_boundary("Bottom",2,0.)
fluid.set_strong_boundary("PtFix",3,0.)
fluid.set_strong_boundary("Top",3,0.)
fluid.solution()[:,3] = rho*g*(fluid.coordinates()[:,1]-0.025)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.export_vtk(outputdir,0,0)
#Computation loop
outf = 50
outf = 1
t = 0.
ii = 0
......@@ -104,7 +101,8 @@ while t < tEnd :
forces = fluid.compute_node_force(dt)
#Computation of the new velocities
vn = p.velocity() + forces * dt / p.mass()
vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
vmax = np.max(np.linalg.norm(vn))
#number of sub time step
nsub = max(1, int(np.ceil((vmax * dt * 4)/min(p.r()))))
print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
......
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