Commit 8311d4ca authored by Nathan Coppin's avatar Nathan Coppin
Browse files

adding doc+adaptingtestcases

parent dd9b6cae
Pipeline #7297 passed with stage
in 2 minutes and 43 seconds
......@@ -223,30 +223,29 @@ class ParticleProblem :
self._lib.particleProblemGetTagId.restype = c_size_t
tag = self._lib.particleProblemGetTagId(self._p, tag.encode())
disks = self.disks()
disks[disks["tag"]==tag]["v"][:] = v
disks["v"][disks["tag"]==tag] = v
segments = self.segments()
segments[segments["tag"]==tag]["v"][:] = v
segments["v"][segments["tag"]==tag] = v
if self._dim == 3 :
triangles = self.triangles()
triangles[triangles["tag"]==tag]["v"][:] = v
triangles["v"][triangles["tag"]==tag] = v
def set_tangent_boundary_velocity(self, tag, vt):
""" Sets the tangent velocity of the frictional boundary.
Only in 2D
Keyword arguments:
tag -- Name of the boundary
vt -- velocity in the tangent direction
"""
assert self._friction_enabled
assert self._dim==2
self._lib.particleProblemGetTagId.restype = c_size_t
tag = self._lib.particleProblemGetTagId(self._p, tag.encode())
disks = self.disks()
segments = self.segments()
disks[disks["tag"]==tag]["vt"] = vt
segments[segments["tag"]==tag]["vt"] = vt
if self._dim == 3 :
triangles = self.triangles()
triangles[triangles["tag"]==tag]["vt"][:] = vt
disks["vt"][disks["tag"]==tag] = vt
segments["vt"][segments["tag"]==tag] = vt
def set_friction_coefficient(self, mu, mat0="default",mat1="default") :
""" Sets the friction coefficient between two materials.
......
import numpy as np
def _advance_particles(particles, f, dt, min_nsub,contact_tol,iteration=0,after_sub_iter=None) :
def _advance_particles(particles, f, dt, min_nsub,contact_tol,iteration=0,after_sub_iter=None,max_split=1) :
"""Contact solver for the grains
Keyword arguments:
......@@ -10,6 +10,8 @@ def _advance_particles(particles, f, dt, min_nsub,contact_tol,iteration=0,after_
min_nsub -- minimal nsub for the particles iterations
contact_tol -- tolerance for the contact solver
iteration -- internal argument -- used to keep track of the number of time step splitting
after_sub_iter -- callback to execute once a sub iteration has been made
max_split -- maximum number of times the time step can be further split if conergence is not reached
"""
# Compute free velocities of the grains
if f is None:
......@@ -23,7 +25,7 @@ def _advance_particles(particles, f, dt, min_nsub,contact_tol,iteration=0,after_
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 == 1:
if iteration == max_split:
for i in range(nsub) :
particles.iterate(dt/nsub, f, tol=contact_tol,force_motion=1)
if after_sub_iter :
......@@ -40,7 +42,7 @@ def _advance_particles(particles, f, dt, min_nsub,contact_tol,iteration=0,after_
after_sub_iter(dt)
def predictor_corrector_iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e-8, external_particles_forces=None, alpha=0.5,after_sub_iter=None) :
def predictor_corrector_iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e-8, external_particles_forces=None, alpha=0.5,after_sub_iter=None,max_split=1) :
"""Predictor-corrector scheme to solve fluid and grains.
Keyword arguments:
......@@ -51,6 +53,8 @@ def predictor_corrector_iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e
contact_tol -- tolerance for the contact solver
external_particles_forces -- external forces applied on the particles
alpha -- parametre of the predictor-corrector scheme [alpha*f(n)+(1-alpha)*f(n+1)]
after_sub_iter -- callback to execute once a sub iteration has been made
max_split -- maximum number of times the time step can be further split if conergence is not reached
"""
if external_particles_forces is None:
external_particles_forces = np.zeros_like(particles.velocity())
......@@ -67,7 +71,7 @@ def predictor_corrector_iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e
fluid.implicit_euler(dt)
sf1 = np.copy(fluid.solution())
f0 = np.copy(fluid.compute_node_force(dt))+external_particles_forces
_advance_particles(particles,f0,dt,min_nsub,contact_tol)
_advance_particles(particles,f0,dt,min_nsub,contact_tol,max_split=max_split)
# Keep same contact forces and positions while giving to the fluid the predicted solid velocities to compute the correction
particles.position()[:,:] = pp0
particles.contact_forces()[:,:] = cp0
......@@ -88,11 +92,11 @@ def predictor_corrector_iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e
fluid.advance_concentration(dt)
# Reset solid velocities
particles.velocity()[:,:] = vp0
_advance_particles(particles,(alpha*f0+(1-alpha)*f1),dt,min_nsub,contact_tol,after_sub_iter=after_sub_iter)
_advance_particles(particles,(alpha*f0+(1-alpha)*f1),dt,min_nsub,contact_tol,after_sub_iter=after_sub_iter,max_split=max_split)
# Give to the fluid the solid information
fluid.set_particles(particles.mass(), particles.volume(), particles.position(), particles.velocity(), particles.contact_forces())
def iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e-8, external_particles_forces=None, fixed_grains=False, after_sub_iter=None) :
def iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e-8, external_particles_forces=None, fixed_grains=False, after_sub_iter=None,max_split=1) :
"""Migflow solver: the solver type depends on the fluid and particles given as arguments.
Keyword arguments:
......@@ -103,6 +107,8 @@ def iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e-8, external_particl
contact_tol -- tolerance for the contact solver
external_particles_forces -- vector of external forces applied on the particles
fixed_grains -- boolean variable specifying if the grains are fixed in the fluid
after_sub_iter -- callback to execute once a sub iteration has been made
max_split -- maximum number of times the time step can be further split if conergence is not reached
Raises:
ValueError -- fluid and particles cannot be both None
......@@ -120,13 +126,13 @@ def iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e-8, external_particl
if fluid.n_fluids() == 2 :
fluid.advance_concentration(dt)
if particles is not None and fluid is None:
_advance_particles(particles,external_particles_forces,dt,min_nsub,contact_tol,after_sub_iter=after_sub_iter)
_advance_particles(particles,external_particles_forces,dt,min_nsub,contact_tol,after_sub_iter=after_sub_iter,max_split=max_split)
if fluid is not None and particles is not None and particles.r() is not None:
if fixed_grains:
fluid.implicit_euler(dt)
fluid.set_particles(particles.mass(), particles.volume(), particles.position(), particles.velocity(),-fluid.compute_node_force(dt))
else:
predictor_corrector_iterate(fluid,particles,dt,min_nsub,contact_tol,external_particles_forces,after_sub_iter=after_sub_iter)
predictor_corrector_iterate(fluid,particles,dt,min_nsub,contact_tol,external_particles_forces,after_sub_iter=after_sub_iter,max_split=max_split)
def particle_changed(fluid,p) :
co = fluid.old_porosity().copy()
......
......@@ -18,32 +18,50 @@
# You should have received a copy of the GNU Lesser General Public License
# along with this program (see COPYING and COPYING.LESSER files). If not,
# see <http://www.gnu.org/licenses/>.
#!/usr/bin/env python
# TESTCASE DESCRIPTION
# Study of the drag on a fixed sphere in a immersed granular flow
# Study of the drag on a fixed cylinder in a granular flow
from migflow import scontact
from migflow import fluid
from migflow import time_integration
import numpy as np
import os
import time
import shutil
import random
numdep = "007"
num = "05"
def genInitialPosition(filename, r, rhop) :
"""Set all the particles centre positions and create the particles objects to add in the computing structure.
Keyword arguments:
filename -- name of the output file
r -- mean radius of the particles
ly -- domain height
lx -- domain width
rhop -- particles density
"""
#Particles structure builder
p = scontact.ParticleProblem(2, friction_enabled = True,rotation_enabled=True)
#Load mesh.msh file specifying physical boundaries names
p.load_msh_boundaries("mesh" + num + ".msh", ["Inner"], material = "Steel")
p.load_msh_boundaries("mesh" + num + ".msh", ["Right","Left","RightUp","RightDown","LeftUp","LeftDown"], material = "Plexi")
p.load_msh_boundaries("mesh" + num + ".msh", ["LockDown"], material = "Steel")
#Particles centre are placed on a regular grid
# Add a grain at each centre position
for y in np.arange(1+r, 2, 2*r*1.05):
for x in np.arange(-0.295+r*1.05, 0.295-r*1.05, 2*r*1.05):
r1 = np.random.uniform(.95*r, 1.05*r)
p.add_particle((x, y), r1, r1**2 * np.pi * rhop,"Sand");
p.write_vtk(outputdir, 0, 0)
# Physical parameters
rapport = 2
g = np.array([0,-9.81])
r = 3e-3
rhop = 2500
# Numerical parameters
dt = 1e-3 #time step
t = 0 #initial time
tEnd = 20 #final time
tEnd = 35 #final time
ii = 0 #iteration number
outf = 10 #iterations between data frames
outf = 100 #iterations between data frames
# Define output directory
inputdir = "depot/depot" + numdep
outputdir = "dry-" + numdep + "-" + str(rapport).replace(".","_")
outputdir = "test" + num
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
#Creating file for drag data
......@@ -56,52 +74,51 @@ if os.path.exists(filename):
#
# PARTICLE PROBLEM
#
p = scontact.ParticleProblem(2,True, True)
p.read_vtk(inputdir,10)
p2 = scontact.ParticleProblem(2,True, True)
p2.read_vtk(inputdir,10)
p2.remove_particles_flag((p2.position()[:,1] + p2.r()[:,0] <1.2))
p2.remove_particles_flag((p2.position()[:,1] - p2.r()[:,0] >1.15))
p2.position()[:,1] += 0.05
flipped = np.copy(p2.position())
flipped[:,0] = -flipped[:,0]
genInitialPosition(outputdir, r, rhop)
p = scontact.ParticleProblem(2,friction_enabled = True,rotation_enabled=True)
p2 = scontact.ParticleProblem(2,friction_enabled = True,rotation_enabled=True)
p.read_vtk(outputdir, 0)
p2.read_vtk(outputdir, 0)
p2.remove_particles_flag((p2.position()[:,1] + p2.r()[:,0] <1.25))
p2.remove_particles_flag((p2.position()[:,1] + p2.r()[:,0] >1.1))
p2.remove_particles_flag((p2.position()[:,0] + p2.r()[:,0] <0.18))
p2.remove_particles_flag((p2.position()[:,0] + p2.r()[:,0] >-0.18))
p2.position()[:,1] -= 0.5
p.clear_boundaries()
p.load_msh_boundaries(inputdir+"/mesh.msh", ["Inner"],material="Steel")
p.load_msh_boundaries(inputdir+"/mesh.msh", ["RightUp","Right","RightDown","LeftDown","Left","LeftUp"],material="Plexi")
p.load_msh_boundaries(inputdir+"/mesh.msh", ["LockDown","LockUp"],material="Steel")
p.remove_particles_flag((p.position()[:,1] + p.r()[:,0] <1.2))
p.add_particles(p2.position(),p2.r(),p2.mass(),v=p2.velocity(),tag="Sand",forced=p2.forced_flag(),contact_forces=p2.contact_forces())
p.add_particles(p2.position()[:,:]+[0,0.05],p2.r(),p2.mass(),v=p2.velocity(),tag="Sand",forced=p2.forced_flag(),contact_forces=p2.contact_forces())
p.add_particles(p2.position()[:,:]+[0,0.1],p2.r(),p2.mass(),v=p2.velocity(),tag="Sand",forced=p2.forced_flag(),contact_forces=p2.contact_forces())
p.add_particles(p2.position()[:,:]+[0,0.15],p2.r(),p2.mass(),v=p2.velocity(),tag="Sand",forced=p2.forced_flag(),contact_forces=p2.contact_forces())
p.add_particles(p2.position()[:,:]+[0,0.2],p2.r(),p2.mass(),v=p2.velocity(),tag="Sand",forced=p2.forced_flag(),contact_forces=p2.contact_forces())
p.set_friction_coefficient(0.4,"Sand","Sand")
p.set_friction_coefficient(0.3,"Sand","Steel")
p.load_msh_boundaries("mesh" + num + ".msh", ["Inner"], material = "Steel")
p.load_msh_boundaries("mesh" + num + ".msh", ["Right","Left","RightUp","RightDown","LeftUp","LeftDown"], material = "Plexi")
p.load_msh_boundaries("mesh" + num + ".msh", ["LockDown"], material = "Steel")
p.set_friction_coefficient(0.15,"Sand","Sand")
p.set_friction_coefficient(0.2,"Sand","Steel")
p.set_friction_coefficient(0.2,"Sand","Plexi")
p.write_vtk(outputdir, 0, 0)
#
# COMPUTATION LOOP
#
while t < tEnd :
opened = 0
F = np.zeros(2)
def accumulate(F):
F+=np.sum(p.get_boundary_impulsion("Inner"),axis=0)/dt
while t < tEnd :
if t>=5 and opened == 0:
opened = 1
p.clear_boundaries()
p.load_msh_boundaries("mesh" + num + ".msh", ["Inner"],material="Steel")
p.load_msh_boundaries("mesh" + num + ".msh", ["Right","Left","RightUp","RightDown","LeftUp","LeftDown"],material="Plexi")
#Adding new particles
if (max(p.position()[:,1])+max(p.r()) <1.4):
p.add_particles(p2.position()[:,:]+[0,0.2] if random.choice([True,False]) else flipped[:,:]+[0,0.2],p2.r(),p2.mass(),v=p2.velocity(),tag="Sand",forced=p2.forced_flag(),contact_forces=p2.contact_forces())
#Updating data structures
G = np.zeros_like(p.velocity())
G[:,1] = p.mass()[:,0]*g[1]
if (max(p.position()[:,1])+max(p.r()) <0.6) and t > 5 and t < 25:
p.add_particles(p2.position(),p2.r(),p2.mass(),v=p2.velocity(),tag="Sand",forced=p2.forced_flag(),contact_forces=p2.contact_forces())
#Solving the contacts
time_integration.iterate(None,p,dt,min_nsub=5,contact_tol=1e-7,external_particles_forces=G,momentum=momentum)
time_integration.iterate(None,p,dt,min_nsub=1,contact_tol=5e-7,external_particles_forces=g*p.mass(),after_sub_iter= lambda subdt : accumulate(F))
#Removing particles outside the hopper
if t>=5 :
p.remove_particles_flag( (p.position()[:,1] >-0.59))
p.remove_particles_flag( (p.position()[:,1] >-0.6))
#Computing mean velocity and writing to file
speedAbove = np.mean(p.velocity()[(p.position()[:,1] < 0.45)*(p.position()[:,1] > 0.3),1])
speedBelow = np.mean(p.velocity()[(p.position()[:,1] > -0.45)*(p.position()[:,1] < -0.3),1])
with open(filename,"a") as file1:
F = np.sum(p.get_boundary_forces(dt,"Inner"),axis=0)
file1.write(str(F[0])+";"+str(F[1])/+";"
+str(t)+";"+str(0 if np.isnan(speedAbove) else speedAbove)
+";"+str(0 if np.isnan(speedBelow) else speedBelow)+"\n")
file1.write(str(F[0])+";"+str(F[1])+";"
+str(t)+";"+str(0 if np.isnan(speedAbove) else speedAbove)+";"
+str(0 if np.isnan(speedBelow) else speedBelow)+"\n")
t += dt
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
......
......@@ -107,16 +107,16 @@ fluid.export_vtk(outputdir, 0.,0)
#
# COMPUTATION LOOP
#
F = np.zeros(2)
def accumulate(F):
F+=np.sum(p.get_boundary_impulsion("Inner"),axis=0)/dt
while t < tEnd :
#Adding new particles
if (max(p.position()[:,1])+max(p.r()) <1.4):
p.add_particles(p2.position()[:,:]+[0,0.2] if random.choice([True,False]) else flipped[:,:]+[0,0.2],p2.r(),p2.mass(),v=p2.velocity(),tag="Sand",forced=p2.forced_flag(),contact_forces=p2.contact_forces())
#Updating data structures
G = np.zeros_like(p.velocity())
G[:,1] = p.mass()[:,0]*g[1]
#Solving the contacts
time_integration.particle_changed(fluid,p)
time_integration.iterate(fluid,p,dt,min_nsub=5,contact_tol=2.5e-7,external_particles_forces=G)
time_integration.iterate(fluid,p,dt,min_nsub=5,contact_tol=2.5e-7,external_particles_forces=g*p.mass(),after_sub_iter = lambda subdt : accumulate(F))
#Removing particles outside the hopper
p.remove_particles_flag( (p.position()[:,1] >-0.59))
#Computing mean velocity and writing to file
......@@ -124,7 +124,7 @@ while t < tEnd :
speedBelow = np.mean(p.velocity()[(p.position()[:,1] > -0.45)*(p.position()[:,1] < -0.3),1])
#print("Speed Above",speedAbove,"Speed Below",speedBelow)
with open(filename,"a") as file1:
F = np.sum(p.get_boundary_forces(dt,"Inner"),axis=0)+ fluid.boundary_force()[0,:];
F += fluid.boundary_force()[0,:];
file1.write(str(F[0])+";"+str(F[1])+";"+str(t)+";"
+str(0 if np.isnan(speedAbove) else speedAbove)+";"+str(0 if np.isnan(speedBelow) else speedBelow)+"\n")
t += dt
......
......@@ -22,18 +22,15 @@ Physical Line("Inner") = {10, 11, 12, 13};
/* Outer Boundary */
Point(16) = {-0.5,1.5,0,lcout};
Point(17) = {-0.5,0.65,0,lcout};
Point(17) = {-0.5,1.65,0,lcout};
Point(18) = {-lout,wout,0,lcout};
Point(19) = {-lout,-wout,0,lcout};
Point(20) = {-0.045,-0.6,0,lcout};
Point(21) = {0.045,-0.6,0,lcout};
Point(22) = {lout,-wout,0,lcout};
Point(23) = {lout,wout,0,lcout};
Point(24) = {0.5,0.65,0,lcout};
Point(25) = {0.5,1.5,0,lcout};
Point(24) = {0.5,1.65,0,lcout};
Line(16) = {16,17};
Line(17) = {17,18};
Line(18) = {18,19};
Line(19) = {19,20};
......@@ -41,17 +38,16 @@ Line(20) = {20,21};
Line(21) = {21,22};
Line(22) = {22,23};
Line(23) = {23,24};
Line(24) = {24,25};
Line(25) = {25,16};
Line(24) = {24,17};
Physical Line("Right") = {18};
Physical Line("Left") = {22};
Physical Line("RightUp") = {16,17};
Physical Line("RightUp") = {17};
Physical Line("RightDown") = {19};
Physical Line("LeftUp") = {23,24};
Physical Line("LeftUp") = {23};
Physical Line("LeftDown") = {21};
Physical Line("LockDown") = {20};
Physical Line("LockUp") = {25};
Physical Line("LockUp") = {24};
Line Loop(20) = {16, 17, 18, 19, 20, 21, 22, 23, 24, 25};
/* Plane Surface */
......
......@@ -20,13 +20,14 @@
# see <http://www.gnu.org/licenses/>.
#!/usr/bin/env python
# TESTCASE DESCRIPTION
# Flow inside a 3D hopper
# Study of the drag on a fixed cylinder in a granular flow
from migflow import scontact
from migflow import fluid
from migflow import time_integration
import numpy as np
import os
import time
def genInitialPosition(filename, r, ly, lx, rhop) :
num = "05"
def genInitialPosition(filename, r, rhop) :
"""Set all the particles centre positions and create the particles objects to add in the computing structure.
Keyword arguments:
......@@ -39,79 +40,89 @@ def genInitialPosition(filename, r, ly, lx, rhop) :
#Particles structure builder
p = scontact.ParticleProblem(3, friction_enabled = True,rotation_enabled=True)
#Load mesh.msh file specifying physical boundaries names
p.load_msh_boundaries("mesh.msh", ["Inner"], material = "Steel")
p.load_msh_boundaries("mesh.msh", ["Left", "Right","Front", "Rear"], material = "Plexi")
p.load_msh_boundaries("mesh.msh", ["LockDown"], material = "Steel")
p.load_msh_boundaries("mesh" + num + ".msh", ["Inner"], material = "Steel")
p.load_msh_boundaries("mesh" + num + ".msh", ["Right","Left","Front","Rear"], material = "Plexi")
p.load_msh_boundaries("mesh" + num + ".msh", ["LockDown"], material = "Steel")
#Particles centre are placed on a regular grid
# Add a grain at each centre position
for y in np.arange(0.86+r, 2.4, 2*r*1.2):
for z in np.arange(0.+r*1.2, 0.05-r*1.2, 2*r*1.2):
for x in np.arange(-0.5+r*1.2, 0.5-r*1.2, 2*r*1.2):
r1 = np.random.uniform(.95*r, 1.05*r)
p.add_particle((x, y, z), r1, 4./3.* r1**3 * np.pi * rhop,"Sand");
#p.position()[:,1] += 0.05
for y in np.arange(1+r, 2, 2*r*1.05):
for x in np.arange(-0.295+r*1.05, 0.295-r*1.05, 2*r*1.05):
for z in np.arange(0+r*1.05,0.05-1.05*r,2*r*1.05):
r1 = np.random.uniform(.95*r, 1.05*r)
p.add_particle((x, y,z), r1, (4/3)*r1**3 * np.pi * rhop,"Sand");
p.write_vtk(outputdir, 0, 0)
# Physical parameters
g = np.array([0,-9.81,0])
r = 3e-3
rhop = 2500
# Numerical parameters
dt = 1e-3 #time step
t = 0 #initial time
tEnd = 35 #final time
ii = 0 #iteration number
outf = 100 #iterations between data frames
# Define output directory
outputdir = "output"
outputdir = "test" + num
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
#Creating file for drag data
filename = outputdir + "/Drag.csv"
if os.path.exists(filename):
print("Remove old results? (y/n)\n")
decision = input()
if decision.lower() == 'y':
os.remove(filename)
# Physical parameters
g = np.array([0,-9.81,0]) # gravity
rhop = 2500 #particles density
r = 7e-3 #radius
ly = 0.5 #particles area height
lx = 0.5 #particles area width
tEnd = 25 #final time
t = 0
# Numerical parameters
dt = 1e-3 #time step
outf = 10 #number of iterations per output
ii = 0
#
# PARTICLE PROBLEM
#
genInitialPosition(outputdir, r, ly, lx, rhop)
genInitialPosition(outputdir, r, rhop)
p = scontact.ParticleProblem(3,friction_enabled = True,rotation_enabled=True)
p2 = scontact.ParticleProblem(3,friction_enabled = True,rotation_enabled=True)
p.read_vtk(outputdir, 0)
p2.read_vtk(outputdir, 0)
p2.remove_particles_flag((p.position()[:,1] + p.r()[:,0] <1.4))
p2.remove_particles_flag((p2.position()[:,1] + p2.r()[:,0] <1.25))
p2.remove_particles_flag((p2.position()[:,1] + p2.r()[:,0] >1.1))
p2.remove_particles_flag((p2.position()[:,0] + p2.r()[:,0] <0.18))
p2.remove_particles_flag((p2.position()[:,0] + p2.r()[:,0] >-0.18))
p2.position()[:,1] -= 0.5
p.clear_boundaries()
p.load_msh_boundaries("mesh.msh", ["Inner"], material = "Steel")
p.load_msh_boundaries("mesh.msh", ["Left", "Right","Front","Rear"], material = "Plexi")
p.load_msh_boundaries("mesh.msh", ["LockDown"], material = "Steel")
p.set_friction_coefficient(0.3,"Sand","Sand")
p.set_friction_coefficient(0.3,"Sand","Steel")
p.set_friction_coefficient(0.15,"Sand","Plexi")
p.load_msh_boundaries("mesh" + num + ".msh", ["Inner"], material = "Steel")
p.load_msh_boundaries("mesh" + num + ".msh", ["Right","Left","RightUp","RightDown","LeftUp","LeftDown","Front","Rear"], material = "Plexi")
p.load_msh_boundaries("mesh" + num + ".msh", ["LockDown"], material = "Steel")
p.set_friction_coefficient(0.15,"Sand","Sand")
p.set_friction_coefficient(0.2,"Sand","Steel")
p.set_friction_coefficient(0.2,"Sand","Plexi")
p.write_vtk(outputdir, 0, 0)
#
# COMPUTATION LOOP
#
opened = 0
F = np.zeros(3)
def accumulate(F):
F+=np.sum(p.get_boundary_impulsion("Inner"),axis=0)/dt
while t < tEnd :
if t>=5 and opened == 0:
opened = 1
p.clear_boundaries()
p.load_msh_boundaries("mesh.msh", ["Inner"],material="Steel")
p.load_msh_boundaries("mesh.msh", ["Left", "Right","Front","Rear"],material="Plexi")
p.remove_particles_flag( (p.position()[:,1] + p.r()[:,0] >-0.61))
if (max(p.position()[:,1])+max(p.r()) <0.85) and t<=20 and t > 5:
p.load_msh_boundaries("mesh" + num + ".msh", ["Inner"],material="Steel")
p.load_msh_boundaries("mesh" + num + ".msh", ["Right","Left","RightUp","RightDown","LeftUp","LeftDown","Front","Rear"],material="Plexi")
#Adding new particles
if (max(p.position()[:,1])+max(p.r()) <0.6) and t > 5 and t < 25:
p.add_particles(p2.position(),p2.r(),p2.mass(),v=p2.velocity(),tag="Sand",forced=p2.forced_flag(),contact_forces=p2.contact_forces())
G = np.zeros_like(p.velocity())
G[:,1] = g[1]*p.mass()[:,0]
time_integration.iterate(None,p,dt,min_nsub=10,contact_tol=5e-7,external_particles_forces=G)
#Solving the contacts
time_integration.iterate(None,p,dt,min_nsub=1,contact_tol=5e-7,external_particles_forces=g*p.mass(),after_sub_iter= lambda subdt : accumulate(F))
#Removing particles outside the hopper
p.remove_particles_flag( (p.position()[:,1] >-0.6))
#Computing mean velocity and writing to file
speedAbove = np.mean(p.velocity()[(p.position()[:,1] < 0.45)*(p.position()[:,1] > 0.3),1])
speedBelow = np.mean(p.velocity()[(p.position()[:,1] > -0.45)*(p.position()[:,1] < -0.3),1])
with open(filename,"a") as file1:
F = np.sum(p_get_boundary_forces(dt,"Inner"),axis=0)
file1.write(str(F[0])+";"+str(F[1])+";"+str(F[2])+";"+str(t)+";"+
str(np.mean(p.velocity()[(p.position()[:,1] < 0.45)*(p.position()[:,1]>0.25),1]))+"\n")
file1.write(str(F[0])+";"+str(F[1])+";"+str(F[2])+";"
+str(t)+";"+str(0 if np.isnan(speedAbove) else speedAbove)+";"
+str(0 if np.isnan(speedBelow) else speedBelow)+"\n")
t += dt
# Output files writing
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
p.write_vtk(outputdir, ioutput, t)
p.write_vtk(outputdir, ioutput, t)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
print("%i : %.2g/%.2g" % (ii, t, tEnd))
......@@ -84,17 +84,17 @@ fluid.set_open_boundary("LockDown",pressure=rho*g[1]*2.1)
fluid.set_open_boundary("LockUp",pressure=0)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces(),init=True)
fluid.export_vtk(outputdir,0,0)
F= np.zeros(3)
def accumulate(F):
F += np.sum(p.get_boundary_impulsion("Inner"))/dt
while t > tEnd :
p.remove_particles_flag( (p.position()[:,1] + p.r()[:,0] >-0.585))
if (max(p.position()[:,1]+p.r()[:,0]) <=0.97) and t<=20:
p.add_particles(p2.position(),p2.r(),p2.mass(),v=p2.velocity(),tag="Sand",forced=p2.forced_flag(),contact_forces=p2.contact_forces())
G = np.zeros_like(p.velocity())
G[:,1] = g[1]*p.mass()[:,0]
p.add_particles(p2.position(),p2.r(),p2.mass(),v=p2.velocity(),tag="Sand",forced=p2.forced_flag(),contact_forces=p2.contact_forces())
time_integration.particle_changed(fluid,p)
time_integration.iterate(fluid,p,dt,min_nsub=8,contact_tol=1e-5,external_particles_forces=G)
time_integration.iterate(fluid,p,dt,min_nsub=8,contact_tol=1e-5,external_particles_forces=g*p.mass(),after_sub_iter=lambda subdt : accumulate(F))
with open(filename,"a") as file1:
F = np.sum(p.get_boundary_forces(dt,"Inner")) + fluid.boundary_force()[0,:]
F += fluid.boundary_force()[0,:]
file1.write(str(F[0])+";"+str(F[1])+";"+str(F[2])+";"+str(t)+";"+
str(np.mean(p.velocity()[(p.position()[:,1] < 0.45)*(p.position()[:,1]>0.25),1]))
+";"+str(np.mean(p.velocity()[(p.position()[:,1] > -0.45)*(p.position()[:,1]<-0.25),1]))+"\n")
......
......@@ -21,18 +21,15 @@ Line Loop(10) = {10, 11, 12, 13};
/* Outer Boundary */
Point(16) = {-0.5,1.5,0,lcout};
Point(17) = {-0.5,0.65,0,lcout};
Point(17) = {-0.5,1.65,0,lcout};
Point(18) = {-lout,wout,0,lcout};
Point(19) = {-lout,-wout,0,lcout};
Point(20) = {-0.05,-0.6,0,lcout};
Point(21) = {0.05,-0.6,0,lcout};
Point(20) = {-0.025,-0.6,0,lcout};
Point(21) = {0.025,-0.6,0,lcout};
Point(22) = {lout,-wout,0,lcout};
Point(23) = {lout,wout,0,lcout};
Point(24) = {0.5,0.65,0,lcout};
Point(25) = {0.5,1.5,0,lcout};
Point(24) = {0.5,1.65,0,lcout};
Line(16) = {16,17};
Line(17) = {17,18};
Line(18) = {18,19};
Line(19) = {19,20};
......@@ -40,22 +37,21 @@ Line(20) = {20,21};
Line(21) = {21,22};
Line(22) = {22,23};
Line(23) = {23,24};
Line(24) = {24,25};
Line(25) = {25,16};
Line(24) = {24,17};
Line Loop(20) = {16, 17, 18, 19, 20, 21, 22, 23, 24, 25};
Line Loop(20) = {17, 18, 19, 20, 21, 22, 23, 24};