Skip to content
GitLab
Commit cb1ff3d2 authored by Nathan Coppin's avatar Nathan Coppin
Browse files

Better drag.py

parent f21448e9
Pipeline #7042 passed with stage
in 2 minutes and 19 seconds
Hide whitespace changes
Inline Side-by-side
python/time_integration.py
View file @ cb1ff3d2
......@@ -21,10 +21,11 @@ def _advance_particles(particles, f, dt, min_nsub,contact_tol,momentum=None,iter
else:
vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
# Estimation of the solid sub-time steps
print("VMAX",vmax)
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 == 3:
if iteration == 1:
for i in range(nsub) :
particles.iterate(dt/nsub, f, tol=contact_tol,force_motion=1)
return
......
scontact/scontact.c
View file @ cb1ff3d2
......@@ -714,6 +714,7 @@ if (fabs(vt)> ((p->a+1)/p->a)*mu*dp) vt *= ((p->a+1)/p->a)*mu*dp/fabs(vt);
#if DIMENSION==2
coordAdd(&p->velocity[c->o1 * DIMENSION], -ct*(p->a/(p->a+1))*c->a1, t);
p->omega[c->o1] -= (1/p->particles[c->o1].r)*(1/((p->a+1)))*ct*c->a1;
#else
cs -= c->cs;
coordAdd(&p->velocity[c->o1 * DIMENSION], -ct*(p->a/(p->a+1))*c->a1, t);
......
testcases/drag-2d/drag.py
View file @ cb1ff3d2
......@@ -33,16 +33,17 @@ import shutil
import random
# Physical parameters
vit = -0.53 #stream velocity
vit = -0.48 #stream velocity
rhop = 2500 #particles density
rho = 500 #fluid density
nu = 1e-8 #fluid kinematic viscosity
g = np.array([0,-21.81])
# Numerical parameters
dt = 1e-3 #time step
t = 0 #initial time
tEnd = 12 #final time
tEnd = 13 #final time
ii = 0 #iteration number
outf = 10 #iterations between data frames
outf = 1 #iterations between data frames
# Define output directory
outputdir = "ellipseForced"
......@@ -61,23 +62,17 @@ if os.path.exists(filename):
#Injected particles
p = scontact.ParticleProblem(2,True,True)
p.read_vtk("depot/depotEllipse",450)
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.position()[:,0] += 0.53005 -0.15
p.position()[:,1] += 0.646 #+ 2*max(p.r())
p.position()[:,1] -=5e-3
p.remove_particles_flag( (p.position()[:,1] > 0.501))
p.velocity()[:,1] = vit
p.velocity()[:,0] = 0
print(min(p.position()[:,1]))
p.read_vtk("depot/depotEllipse",450)
p.remove_particles_flag( (p.position()[:,1] < -0.2339))
p.position()[:,:] = p.position()[:,(1,0)]
p.position()[:,:] += [0.38005,0.647]
p.remove_particles_flag( (p.position()[:,1] > 0.507))
p.velocity()[:,:] = [0,0]
flipped = np.copy(p.position())
flipped[:,0] = -flipped[:,0]
p.contact_forces()[:,:] = 0
p.omega()[:] = 0
p.forced_flag()[:] = 1
#p.r()[:,0] = 0.999
p.forced_flag()[:] = 1
p.mass()[:] *= 0.05
p.write_vtk(outputdir + "/In_particles", 0, 0)
#Injected particles
......@@ -91,21 +86,24 @@ 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.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] -= 1.165
p2.add_particles(p.position(),p.r(),p.mass(),v=np.zeros_like(p.velocity()),tag="Sand",forced=np.zeros_like(p.forced_flag()),contact_forces=p.contact_forces())
p2.position()[p2.position()[:,1]>0,1] -= 0.88
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.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.add_particles(p.position()[:,:] - [0,0.291],p.r(),p.mass(),v=np.zeros_like(p.velocity()),tag="Sand",forced=np.zeros_like(p.forced_flag()),contact_forces=p.contact_forces())
p2.add_particles(flipped - [0,0.291*2],p.r(),p.mass(),v=np.zeros_like(p.velocity()),tag="Sand",forced=np.zeros_like(p.forced_flag()),contact_forces=p.contact_forces())
p2.add_particles(p.position()[:,:] - [0,0.291*3],p.r(),p.mass(),v=np.zeros_like(p.velocity()),tag="Sand",forced=np.zeros_like(p.forced_flag()),contact_forces=p.contact_forces())
p2.add_particles(flipped - [0,0.291*4],p.r(),p.mass(),v=np.zeros_like(p.velocity()),tag="Sand",forced=np.zeros_like(p.forced_flag()),contact_forces=p.contact_forces())
p2.remove_particles_flag( (p2.position()[:,1]**2 + p2.position()[:,0]**2 > (0.01905 + max(p2.r()))**2))
p2.r()[:] *= 0.99
p2.set_friction_coefficient(0.28,"Sand","Sand")
p2.set_friction_coefficient(0.28,"Sand","Steel")
p2.set_friction_coefficient(0.14,"Sand","Plexi")
p2.set_tangent_boundary_velocity("Left",0)#vit)
p2.set_tangent_boundary_velocity("Right",0)#-vit)
p2.forced_flag()[p2.position()[:,1] + 2*p2.r()[:,0] >= 0.5] = 1
p2.contact_forces()[:,:] = 0
p2.velocity()[:,0] = 0
p2.velocity()[p2.position()[:,1] < 0,1] = 0
p2.omega()[:] = 0
#p2.r()[:,0] *= 0.999
p2.mass()[:,0] *= 0.05
p2.set_get_momentum(1)
p2.write_vtk(outputdir+"/Flow", 0, 0)
......@@ -128,23 +126,36 @@ p2.write_vtk(outputdir+"/Flow", 0, 0)
#
# COMPUTATION LOOP
#++
print(p.n_particles())
print(p2.n_particles())
tic = time.time()
vitC = -0.01
while t < tEnd :
momentum = np.zeros((3,2))
momentum = np.zeros_like(p2.momentum())
if vitC >= vit and t>=500*dt:
vitC += 2*dt*vit
vitB = vitC#*0.72 if int((t+dt)/dt)%2==0 else vitC#-0.002 if int((t/dt))%3==0 else vitC+0.002
vitC2 = vitC#*72.0 if int(t/dt)%2==0 else vitC
print("vitC", vitC, "vitB", vitB)
p2.forced_flag()[p2.position()[:,1] + 2*p2.r()[:,0] < 0.5] = 0
p2.forced_flag()[(p2.position()[:,1] + 2*p2.r()[:,0] <= -0.5)+(p2.position()[:,1] + 2*p2.r()[:,0] >= 0.5)] = 1
p2.velocity()[p2.forced_flag()==1,:] = [0, vitC2]
p2.velocity()[(p2.position()[:,1] + 2*p2.r()[:,0] <= -0.5),1] = vitB
p2.omega()[p2.forced_flag()==1] = 0
if max(p2.position()[:,1]) < 0.5 :
p2.add_particles(p.position(),p.r(),p.mass(),v=p.velocity(),tag="Sand",forced=p.forced_flag(),contact_forces=p.contact_forces())
p2.add_particles(p.position() if random.choice([True,False]) else flipped,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,10 if t <= 100*dt else 10,contact_tol=1.5e-8,momentum=momentum)
G = np.zeros_like(p2.velocity())
G[:,0] = p2.mass()[:,0]*g[1]
time_integration.iterate(None,p2,dt,20,contact_tol=3e-8,momentum=momentum)
#momentum += fluid.boundary_force()[0:3,:]
print(momentum)
p2.forced_flag()[p2.position()[:,1] + 2*p2.r()[:,0] < 0.5] = 0
p2.forced_flag()[p2.position()[:,1] + 2*p2.r()[:,0] < -0.5] = 1
p2.velocity()[p2.position()[:,1] + 2*p2.r()[:,0] < -0.5,:] = [0,vit]
p2.omega()[p2.position()[:,1] + 2*p2.r()[:,0] < -0.5] = 0
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:
file1.write(str(momentum[0,0]) + ";" + str(momentum[0,1]) + ";" + str(t) + ";" + str(vit) + "\n")
file1.write(str(momentum[0,0]) + ";" + str(momentum[0,1]) + ";" + str(momentum[1,0]) + ";"
+ str(momentum[1,1]) + ";" + str(momentum[2,0]) + ";" + str(momentum[2,1]) + ";"
+ 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 :
......
testcases/drag-2d/plot.py
View file @ cb1ff3d2
......@@ -6,48 +6,44 @@ import random
import csv
import matplotlib.pyplot as plt
filename = "ellipseGravity/Drag.csv"
filename = "ellipseForced/Drag.csv"
with open(filename) as file1:
reader = csv.reader(file1,delimiter = ";",quoting=csv.QUOTE_NONNUMERIC)
data = [r for r in reader]
data = np.array(data);
react = data[:,5]#(data[:,0]**2 + data[:,1]**2)**0.5
#react2 = data[:,3]
#react3 = data[:,5]
react = data[:,1]
react2 = data[:,3]
react3 = data[:,5]
prod = np.zeros_like(react)
#prod2 = np.zeros_like(react2)
#prod3 = np.zeros_like(react3)
ranges = 100
prod2 = np.zeros_like(react2)
prod3 = np.zeros_like(react3)
ranges = 500
print(data)
for i in range(len(react)):
if i <ranges:
prod[i] = np.sum(react[0:i+ranges])/(i+ranges)
#prod2[i] = np.sum(react2[0:i+ranges])/(i+ranges)
#prod3[i] = np.sum(react3[0:i+ranges])/(i+ranges)
prod2[i] = np.sum(react2[0:i+ranges])/(i+ranges)
prod3[i] = np.sum(react3[0:i+ranges])/(i+ranges)
elif len(react)-i<ranges:
prod[i] = np.sum(react[i-ranges:])/((len(react)-i)+ranges)
#prod2[i] = np.sum(react2[i-ranges:])/((len(react)-i)+ranges)
#prod3[i] = np.sum(react3[i-ranges:])/((len(react)-i)+ranges)
prod2[i] = np.sum(react2[i-ranges:])/((len(react)-i)+ranges)
prod3[i] = np.sum(react3[i-ranges:])/((len(react)-i)+ranges)
else:
prod[i] = np.sum(react[i-ranges:i+ranges])/(2*ranges)
#prod2[i] = np.sum(react2[i-ranges:i+ranges])/(2*ranges)
#prod3[i] = np.sum(react3[i-ranges:i+ranges])/(2*ranges)
prod2[i] = np.sum(react2[i-ranges:i+ranges])/(2*ranges)
prod3[i] = np.sum(react3[i-ranges:i+ranges])/(2*ranges)
ProdToMean = prod#[5500:15000]
print(np.mean(ProdToMean[2500:]))
mean = np.ones_like(ProdToMean)*np.mean(ProdToMean[5500:12000])
TimeVel = data[2000:14000, 7]
print(np.mean(prod[5000:]))
mean = np.ones_like(prod[5000:])*np.mean(prod[5000:])
TimeVel = (data[5000:, 7] + data[5000:, 8])/2
plt.plot(data[:,6],-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[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.plot(data[:,6],-prod2, color = 'purple')
plt.plot(data[:,6],-prod3, color = 'green')
plt.plot(data[5000:,6],-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[2500:])*1000))
plt.ylim([-1,360])
plt.xlim([.5,24])
plt.title("v = %.1f mm/s"%(np.mean(TimeVel)*1000))
plt.ylim([-1,150])
plt.xlim([.0,20])
plt.show()
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