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

Better managemetn of time step splitting for force recup

parent 856f0d0d
......@@ -37,7 +37,7 @@ def _advance_particles(particles, f, dt, min_nsub,contact_tol,momentum=None,iter
if momentum is not None:
#print(particles.momentum()*nsub/dt)
momentum += particles.momentum()
if momentum is not None:
if momentum is not None and iteration == 0:
momentum /= dt
......
......@@ -33,19 +33,19 @@ import shutil
import random
# Physical parameters
vit = -0.07 #stream velocity
vit = -0.6 #stream velocity
rhop = 2500 #particles density
rho = 500 #fluid density
nu = 1e-8 #fluid kinematic viscosity
# Numerical parameters
dt = 1e-3 #time step
t = 0 #initial time
tEnd = 12.5 #final time
tEnd = 25 #final time
ii = 0 #iteration number
outf = 10 #iterations between data frames
# Define output directory
outputdir = "vit007"
outputdir = "vit007fois2+vitpar"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
filename = outputdir + "/Drag.csv"
......@@ -61,7 +61,7 @@ if os.path.exists(filename):
#Injected particles
p = scontact.ParticleProblem(2,True,True)
p.read_vtk("depot/outputCompactPeriodic2",985)
p.read_vtk("depot/outputCompactPeriodicPoly2",220)
#p.remove_particles_flag( (p.position()[:,1] <= -0.231 - p.r()[:,0]))
#pos = np.zeros_like(p.position())
#pos[:,0] = p.position()[:,1]
......@@ -86,19 +86,19 @@ p.write_vtk(outputdir + "/In_particles", 0, 0)
#Flow particles
p2 = scontact.ParticleProblem(2,True,True)
p2.read_vtk("depot/output",1800)
p2.read_vtk(outputdir + "/In_particles",0)
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.5-2*p2.r()[:,0]))
#p2.remove_particles_flag( (p2.position()[:,1] < 0.5-2*p2.r()[:,0]))
#p2.add_particles(p.position(),p.r(),p.mass(),v=p.velocity(),tag="Sand",forced=p.forced_flag(),contact_forces=p.contact_forces())
#p2.position()[:,1] -= 0.9
p2.position()[:,1] -= 0.9
#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.4,"Sand","Sand")
p2.set_friction_coefficient(0.4,"Sand","Steel")
p2.set_friction_coefficient(0.05,"Sand","Plexi")
#p2.set_tangent_boundary_velocity("Left",-vit)
#p2.set_tangent_boundary_velocity("Right",vit)
p2.set_tangent_boundary_velocity("Left",-vit)
p2.set_tangent_boundary_velocity("Right",vit)
p2.forced_flag()[(p2.position()[:,1] + 2*p2.r()[:,0] < 0.5)*(p2.position()[:,1] + 2*p2.r()[:,0] >= 0.4)] = 1
p2.velocity()[:,1] = vit
p2.contact_forces()[:,:] = 0
......@@ -136,7 +136,7 @@ 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.551))
#time_integration.particle_changed(fluid,p2)
time_integration.iterate(None,p2,dt,min_nsub=40 if t < 1000*dt else 10,contact_tol=1e-8,momentum=momentum)
time_integration.iterate(None,p2,dt,80 if t <= 100*dt else 10,contact_tol=1e-8,momentum=momentum)
#momentum += fluid.boundary_force()[0:3,:]
print(momentum)
p2.forced_flag()[p2.position()[:,1] + 2*p2.r()[:,0] < 0.5] = 0
......
......@@ -6,7 +6,7 @@ import random
import csv
import matplotlib.pyplot as plt
filename = "output009/Drag.csv"
filename = "vit007fois2+vitpar/Drag.csv"
with open(filename) as file1:
reader = csv.reader(file1,delimiter = ";",quoting=csv.QUOTE_NONNUMERIC)
data = [r for r in reader]
......@@ -17,7 +17,7 @@ react = data[:,1]#(data[:,0]**2 + data[:,1]**2)**0.5
prod = np.zeros_like(react)
#prod2 = np.zeros_like(react2)
#prod3 = np.zeros_like(react3)
ranges = 500
ranges = 100
print(data)
for i in range(len(react)):
if i <ranges:
......@@ -34,20 +34,20 @@ for i in range(len(react)):
#prod3[i] = np.sum(react3[i-ranges:i+ranges])/(2*ranges)
ProdToMean = prod#[5500:15000]
print(np.mean(ProdToMean[11000:12500]))
mean = np.ones_like(ProdToMean)*np.mean(ProdToMean[6000:15000])
TimeVel = data[:, 3]
plt.plot(data[:,6],-prod,color='black')
print(np.mean(ProdToMean[2000:8000]))
mean = np.ones_like(ProdToMean)*np.mean(ProdToMean[2000:8000])
#TimeVel = data[6000:14000, 2]
plt.plot(data[:,2],-prod,color='black')
#plt.plot(data[:,6],-prod2, color = 'purple')
#plt.plot(data[:,6],-prod3, color = 'green')
plt.plot(data[:,6],-mean,color = 'red')
plt.plot(data[:,2],-mean,color = 'red')
#plt.plot(data[np.argwhere(prod<=-100),6],-mean,color = 'red')
#mean = np.ones_like(data)*np.mean(data[:,1])
#plt.plot(data[:,2],data[:,1])
#plt.plot(data[:,2],mean,color = 'red')
plt.xlabel("Time [s]",size=15)
plt.ylabel("Drag force [N]",size=15)
#plt.title("v = %.1f mm/s")%(np.mean(TimeVel[5500:15500])*1000),size=19)
#plt.title("v = %.1f mm/s"%(np.mean(TimeVel[5500:14500])*1000))
plt.ylim([-1,360])
#plt.xlim([1.5,24])
plt.xlim([.5,24])
plt.show()
......@@ -50,7 +50,7 @@ def genInitialPosition(filename, r, ly, lx, rhop) :
p.load_msh_boundaries("mesh.msh", ["LockDown"], material = "Steel")
#P articles centre are placed on a regular grid
# Add a grain at each centre position
for y in np.arange(0.75, 1.55, 2*r*1.2):
for y in np.arange(0.86+r, 1.60, 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)
......@@ -72,14 +72,14 @@ if os.path.exists(filename):
# Physical parameters
g = np.array([0,-9.81,0]) # gravity
rhop = 2500 #particles density
r = 1e-2 #radius
r = 3e-3 #radius
ly = 0.5 #particles area height
lx = 0.5 #particles area width
tEnd = 30 #final time
tEnd = 25 #final time
t = 0
# Numerical parameters
dt = 1e-3 #time step
outf = 1 #number of iterations per output
outf = 100 #number of iterations per output
ii = 0
#
# PARTICLE PROBLEM
......@@ -96,9 +96,9 @@ p.load_msh_boundaries("mesh.msh", ["Left", "Right","Front","Rear"], material = "
p.load_msh_boundaries("mesh.msh", ["LockDown"], material = "Steel")
p.set_get_momentum(1)
p.set_friction_coefficient(0.13,"Sand","Sand")
p.set_friction_coefficient(0.13,"Sand","Steel")
p.set_friction_coefficient(0.03,"Sand","Plexi")
p.set_friction_coefficient(0.4,"Sand","Sand")
p.set_friction_coefficient(0.4,"Sand","Steel")
p.set_friction_coefficient(0.2,"Sand","Plexi")
print("r = %g, m = %g\n" % (p.r()[0], p.mass()[0]))
print("RHOP = %g" % rhop)
......@@ -106,29 +106,26 @@ print("Number of particles : \n")
print(p.n_particles())
p.write_vtk(outputdir, 0, 0)
G = np.zeros_like(p.velocity())
G[:,1] = g[1]*p.mass()[:,0]
opened = 0
tic = time.time()
while t < tEnd :
if t>=2.5 and opened == 0:
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")
if np.any(p.position()[:,1]+ p.r()[:,0] < -0.61):
p.remove_particles_flag( (p.position()[:,1] + p.r()[:,0] >-0.61))
G = np.zeros_like(p.velocity())
G[:,1] = g[1]*p.mass()[:,0]
if (max(p.position()[:,1])+max(p.r()) <0.85) and t<=15 and t > 2.5:
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.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]
G = np.zeros_like(p.velocity())
G[:,1] = g[1]*p.mass()[:,0]
momentum = np.zeros_like(p.momentum())
time_integration.iterate(None,p,dt,min_nsub=10,contact_tol=2e-8,external_particles_forces=G,momentum=momentum)
time_integration.iterate(None,p,dt,min_nsub=10,contact_tol=5e-9,external_particles_forces=G,momentum=momentum)
with open(filename,"a") as file1:
file1.write(str(momentum[0,0]) + ";" + str(momentum[0,1])+ ";" + str(momentum[0,2]) + ";" + str(t) + ";" + str(np.mean(p.velocity()[(p.position()[:,1] < 0.4)*(p.position()[:,1]>0.2),1])) +"\n")
file1.write(str(momentum[0,0]) + ";" + str(momentum[0,1])+ ";" + str(momentum[0,2]) + ";" + str(t) + ";" + str(np.mean(p.velocity()[(p.position()[:,1] < 0.45)*(p.position()[:,1]>0.25),1])) +"\n")
t += dt
# Output files writing
if ii %outf == 0 :
......
......@@ -123,7 +123,7 @@ G[:,1] = p.mass()[:,0]*g[1]
while t < tEnd :
#print(p.omega())
print(np.sum(p.omega()**2)**0.5)
time_integration.iterate(None, p, dt, min_nsub=15, external_particles_forces=G)
time_integration.iterate(fluid, p, dt, min_nsub=15, external_particles_forces=G)
t += dt
#Output files writting
if ii %outf == 0 :
......
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