Commit d3feed44 authored by Nathan Coppin's avatar Nathan Coppin
Browse files

polydisperse testcase in drag

parent e6fc4ff7
Pipeline #6758 passed with stage
in 2 minutes and 51 seconds
......@@ -24,7 +24,7 @@ def _advance_particles(particles, f, dt, min_nsub,contact_tol,momentum=None,iter
nsub = max(min_nsub, int(np.ceil((vmax * dt * 8)/min(particles.r()))))
print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(particles.r()))
#If time step was split too many times, ignore the convergence and move the grains
if iteration == 10:
if iteration == 5:
for i in range(nsub) :
particles.iterate(dt/nsub, f, tol=contact_tol,force_motion=1)
return
......
......@@ -33,7 +33,7 @@ import shutil
import random
# Physical parameters
vit = -0.05 #stream velocity
vit = -0.15 #stream velocity
rhop = 2500 #particles density
rho = 1000 #fluid density
nu = 1e-6 #fluid kinematic viscosity
......@@ -45,7 +45,7 @@ ii = 0 #iteration number
outf = 1 #iterations between data frames
# Define output directory
outputdir = "again81"
outputdir = "again91"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
filename = outputdir + "/Drag.csv"
......@@ -61,12 +61,13 @@ if os.path.exists(filename):
#Injected particles
p = scontact.ParticleProblem(2,True,True)
p.read_vtk("depot/outputPeriodic",435)
p.read_vtk("depot/outputPolydispersePeriodic",435)
p.remove_particles_flag( (p.position()[:,1] <= -0.231 - p.r()[:,0]))
pos = np.zeros_like(p.position())
pos[:,0] = p.position()[:,1]
pos[:,1] = p.position()[:,0]
p.position()[:,:] = pos[:,:]
p.mass()[:,0] *= 0.05*2/3
p.position()[:,0] += 0.53005 -0.15
p.position()[:,1] += .55 #+ 2*max(p.r())
p.remove_particles_flag( (p.position()[:,1] > 0.405))
......@@ -74,23 +75,30 @@ p.position()[:,1] +=2.5e-3
p.velocity()[:,1] = vit
p.velocity()[:,0] = 0
p.contact_forces()[:,:] = 0
p.forced_flag()[:] = 1
p.omega()[:] = 0
p.forced_flag()[:] = 1
#p.r()[:,0] = 0.999
p.write_vtk(outputdir + "/In_particles", 0, 0)
#Injected particles
#Flow particles
p2 = scontact.ParticleProblem(2,True,True)
p2.read_vtk("depot/output",1800)
p2.read_vtk("depot/outputPolydisperse",1800)
p2.clear_boundaries()
p2.load_msh_boundaries("mesh.msh", ["Inner"],material = "Steel")
p2.load_msh_boundaries("mesh.msh", ["Left", "Right"],material = "Plexi")
p2.remove_particles_flag( (p2.position()[:,1] < 0.4-p2.r()[:,0]))
p2.set_friction_coefficient(0.,"Sand","Sand")
p2.set_friction_coefficient(0.,"Sand","Steel")
p2.set_tangent_boundary_velocity("Left",vit)
p2.set_tangent_boundary_velocity("Right",-vit)
#p2.add_particles(p.position(),p.r(),p.mass(),v=np.zeros_like(p.velocity()),tag="Sand",forced=p.forced_flag(),contact_forces=p.contact_forces())
p2.set_friction_coefficient(0.2,"Sand","Sand")
p2.set_friction_coefficient(0.2,"Sand","Steel")
p2.set_tangent_boundary_velocity("Left",-vit)
p2.set_tangent_boundary_velocity("Right",vit)
p2.contact_forces()[:,:] = 0
p2.velocity()[:,:] = 0
p2.omega()[:] = 0
p2.r()[:,0] *= 0.999
p2.mass()[:,0] *= 0.05*2/3
p2.set_get_momentum(1)
p2.write_vtk(outputdir+"/Flow", 0, 0)
......@@ -112,6 +120,7 @@ p2.write_vtk(outputdir+"/Flow", 0, 0)
#
# COMPUTATION LOOP
#++
print(p.n_particles())
tic = time.time()
while t < tEnd :
momentum = np.zeros((3,2))
......@@ -119,12 +128,13 @@ while t < tEnd :
p2.add_particles(p.position(),p.r(),p.mass(),v=p.velocity(),tag="Sand",forced=p.forced_flag(),contact_forces=p.contact_forces())
p2.remove_particles_flag( (p2.position()[:,1] + p2.r()[:,0] >-0.451))
#time_integration.particle_changed(fluid,p2)
time_integration.iterate(None,p2,dt,min_nsub=5,contact_tol=1e-7,momentum=momentum)
time_integration.iterate(None,p2,dt,min_nsub=40 if t < 1000*dt else 10,contact_tol=1e-8,momentum=momentum)
#momentum += fluid.boundary_force()
print(momentum)
p2.forced_flag()[p2.position()[:,1] + 2*p2.r()[:,0] < 0.4] = 0
p2.forced_flag()[p2.position()[:,1] + 2*p2.r()[:,0] < -0.4] = 1
p2.velocity()[p2.position()[:,1] + 2*p2.r()[:,0] < -0.4,:] = [0,vit]
p2.omega()[p2.position()[:,1] + 2*p2.r()[:,0] < -0.4] = 0
with open(filename,"a") as file1:
file1.write(str(momentum[0,0]) + ";" + str(momentum[0,1]) + ";" + str(t) + ";" + str(vit) + "\n")
t += dt
......
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