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

Corrected index problem in getting mu for disks

parent 459b46e3
Pipeline #6656 passed with stage
in 2 minutes and 49 seconds
......@@ -296,6 +296,7 @@ class ParticleProblem :
def read_vtk(self,dirname,i,contact_factor=1) :
"""Read an existing output file to set the particle data."""
self._lib.particleProblemClearBoundaries(self._p)
x,_,d,cdata,fdata = VTK.read(dirname+"/particles_%05d.vtp"%i)
matnames = VTK.string_array_decode(fdata["MaterialNames"])
matmap = {}
......@@ -311,7 +312,7 @@ class ParticleProblem :
tagmap = {}
for j,n in enumerate(tagnames) :
tagmap[j] = self._lib.particleProblemGetTagId(self._p,n)
self._lib.particleProblemClearBoundaries(self._p)
offsets = np.hstack([[0],cells["offsets"]])
el = cells["connectivity"]
tags = np.vectorize(tagmap.get)(cdata["Tag"])
......
......@@ -32,7 +32,7 @@ def _advance_particles(particles, f, dt, min_nsub,contact_tol,momentum=None,iter
for i in range(nsub) :
if not particles.iterate(dt/nsub, f, tol=contact_tol):
print("Splitting time-step to level %d..."%(iteration+1))
_advance_particles(particles,f,dt/nsub,2,contact_tol/2,iteration+1,momentum=momentum)
_advance_particles(particles,f,dt/nsub,2,contact_tol/2,momentum=momentum,iteration=iteration+1)
print("Convergence reached for subinvervals of level %d"%(iteration+1))
if momentum is not None:
momentum += particles.momentum()*nsub/dt
......
......@@ -448,6 +448,7 @@ static int contactParticleParticleSolve(Contact *c, ParticleProblem *p, double d
const double mu = particleProblemGetMu(p,p->particleMaterial[c->o0],p->particleMaterial[c->o1]);
#if DIMENSION==2
if(isnan(dp)){printf("hello\n");}
if (fabs(vt)> ((p->a+1)/p->a)*mu*dp) vt *= ((p->a+1)/p->a)*mu*dp/fabs(vt);
ct = vt;
if (dp == 0.){
......@@ -566,8 +567,9 @@ static int contactParticleDiskSolve(Contact *c, ParticleProblem *p, double dt, d
#endif
}
#endif
const double mu = particleProblemGetMu(p,p->particleMaterial[c->o1],p->segmentMaterial[c->o0]);
const double mu = particleProblemGetMu(p,p->particleMaterial[c->o1],p->diskMaterial[c->o0]);
#if DIMENSION==2
if(isnan(dp)){printf("hello\n");}
if (fabs(vt)> ((p->a+1)/p->a)*mu*dp) vt *= ((p->a+1)/p->a)*mu*dp/fabs(vt);
ct = vt;
if (dp == 0.){
......@@ -671,6 +673,7 @@ static int contactParticleSegmentSolve(Contact *c, ParticleProblem *p, double dt
#endif
const double mu = particleProblemGetMu(p,p->particleMaterial[c->o1],p->segmentMaterial[c->o0]);
#if DIMENSION==2
if(isnan(dp)){printf("hello\n");}
if (fabs(vt)> ((p->a+1)/p->a)*mu*dp) vt *= ((p->a+1)/p->a)*mu*dp/fabs(vt);
ct = vt;
if (dp == 0.){
......@@ -1266,7 +1269,7 @@ void _get_momentum(ParticleProblem *p) {
#if FRICTION_ENABLED && DIMENSION==2
p->momentum[DIMENSION*p->diskTag[p->contacts[i].o0] +j] += t[j]*p->particles[p->contacts[i].o1].m*p->contacts[i].ct;
#elif FRICTION_ENABLED && DIMENSION==3
p->momentum[DIMENSION*p->diskTag[p->contacts[i].o0] +j] += s[j]*p->particles[p->contacts[i].o1].m*p->contacts[i].ct;
p->momentum[DIMENSION*p->diskTag[p->contacts[i].o0] +j] += t[j]*p->particles[p->contacts[i].o1].m*p->contacts[i].ct;
p->momentum[DIMENSION*p->diskTag[p->contacts[i].o0] +j] += s[j]*p->particles[p->contacts[i].o1].m*p->contacts[i].cs;
#endif
}
......@@ -1425,7 +1428,11 @@ void particleProblemRemoveParticles(ParticleProblem *p, const int *keep_flag) {
vectorRemoveFlag(p->particleMaterial,keep_flag,1);
vectorRemoveFlag(p->ForcedFlag,keep_flag,1);
#if FRICTION_ENABLED && ROTATION_ENABLED
vectorRemoveFlag(p->omega,keep_flag,DIMENSION==2 ? 1 : 3);
#if DIMENSION==2
vectorRemoveFlag(p->omega,keep_flag,1);
#else
vectorRemoveFlag(p->omega,keep_flag,3);
#endif
#endif
// update contacts
size_t * map = malloc(sizeof(size_t)*vectorSize(p->particles));
......
......@@ -33,7 +33,7 @@ import shutil
import random
# Physical parameters
vit = -0.02 #stream velocity
vit = -0.05 #stream velocity
rhop = 2500 #particles density
rho = 1000 #fluid density
nu = 1e-6 #fluid kinematic viscosity
......@@ -42,7 +42,7 @@ dt = 1e-3 #time step
t = 0 #initial time
tEnd = 50 #final time
ii = 0 #iteration number
outf = 50 #iterations between data frames
outf = 1 #iterations between data frames
# Define output directory
outputdir = "again81"
......@@ -60,16 +60,17 @@ if os.path.exists(filename):
#
#Injected particles
p = scontact.ParticleProblem(2)
p.read_vtk("depot/output9",10)
p.remove_particles_flag( (p.position()[:,1] < -p.r()[:,0]))
p = scontact.ParticleProblem(2,True,True)
p.read_vtk("depot/outputPeriodic",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.position()[:,0] += 0.15005
p.position()[:,1] += 0.35
p.remove_particles_flag( (p.position()[:,1] > 0.25+np.max(p.r())))
p.position()[:,0] += 0.53005 -0.15
p.position()[:,1] += .55 #+ 2*max(p.r())
p.remove_particles_flag( (p.position()[:,1] > 0.405))
p.position()[:,1] +=2.5e-3
p.velocity()[:,1] = vit
p.velocity()[:,0] = 0
p.contact_forces()[:,:] = 0
......@@ -79,48 +80,51 @@ p.write_vtk(outputdir + "/In_particles", 0, 0)
#Flow particles
p2 = scontact.ParticleProblem(2,True)
p2.load_msh_boundaries("mesh2.msh", ["Inner"],material = "Steel")
p2.load_msh_boundaries("mesh2.msh", ["Left", "Right","Bottom", "Top"],material = "Plexi")
p2.set_friction_coefficient(0.2,"Sand","Sand")
p2.set_friction_coefficient(0.15,"Sand","Steel")
p2.set_friction_coefficient(0.15,"Sand","Plexi")
p2 = scontact.ParticleProblem(2,True,True)
p2.read_vtk("depot/output",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.contact_forces()[:,:] = 0
p2.set_get_momentum(1)
p2.add_particles(p.position(),p.r(),p.mass(),v=p.velocity(),tag="Sand",forced=p.forced_flag(),contact_forces=p.contact_forces())
p2.write_vtk(outputdir+"/Flow", 0, 0)
#
# FLUID PROBLEM
#
fluid = fluid.FluidProblem(2,0,[nu*rho],[rho])
#fluid = fluid.FluidProblem(2,0,[nu*rho],[rho])
#Set the mesh geometry for the fluid computation
fluid.load_msh("mesh2.msh")
fluid.set_wall_boundary("Inner")
fluid.set_wall_boundary("Left", velocity = [0,vit])
fluid.set_wall_boundary("Right", velocity = [0,vit])
fluid.set_open_boundary("Bottom", velocity = [0,vit])
fluid.set_open_boundary("Top",pressure = 0)
fluid.set_particles(p2.mass(), p2.volume(), p2.position(), p2.velocity(), p2.contact_forces(),init=True)
fluid.export_vtk(outputdir+"/Flow", 0.,0)
#fluid.load_msh("mesh2.msh")
#fluid.set_wall_boundary("Inner")
#fluid.set_wall_boundary("Left", velocity = [0,vit])
#fluid.set_wall_boundary("Right", velocity = [0,vit])
#fluid.set_open_boundary("Bottom", velocity = [0,vit])
#fluid.set_open_boundary("Top",pressure = 0)
#fluid.set_particles(p2.mass(), p2.volume(), p2.position(), p2.velocity(), p2.contact_forces(),init=True)
#fluid.export_vtk(outputdir+"/Flow", 0.,0)
#
# COMPUTATION LOOP
#
#++
tic = time.time()
while t < tEnd :
momentum = np.zeros((5,2))
if max(p2.position()[:,1]) < 0.25 + p2.r()[0,0] :
momentum = np.zeros((3,2))
if max(p2.position()[:,1]) < 0.405 :
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.1))
time_integration.particle_changed(fluid,p2)
time_integration.iterate(fluid,p2,dt,min_nsub=5,contact_tol=5e-9,momentum=momentum)
momentum += fluid.boundary_force()
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)
#momentum += fluid.boundary_force()
print(momentum)
p2.forced_flag()[p2.position()[:,1] + 2*p2.r()[:,0] < 0.2] = 0
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]
with open(filename,"a") as file1:
file1.write(str(momentum[0,0]) + ";" + str(momentum[0,1]) + ";" + str(t) + ";" + str(vit) + "\n")
t += dt
......@@ -128,7 +132,7 @@ while t < tEnd :
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
p2.write_vtk(outputdir+"/Flow", ioutput, t)
fluid.export_vtk(outputdir+"/Flow", t, ioutput)
#fluid.export_vtk(outputdir+"/Flow", t, ioutput)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
rin = 0.0064;
wout = 0.125;
lout1 = 0.02228;
lout2 = 0.03022;
rin = 0.01905;
wout = 1.;
lout = 0.15;
lcin = wout/20;
lcout = wout/20;
......@@ -22,18 +21,24 @@ Line Loop(10) = {10, 11, 12, 13};
Physical Line("Inner") = {10, 11, 12, 13};
/* Outer Boundary */
Point(20) = {wout/5, -wout, 0, lcout};
Point(21) = {wout/5, wout, 0, lcout};
Point(22) = {-wout/5, wout, 0, lcout};
Point(23) = {-wout/5, -wout, 0, lcout};
Point(20) = {lout, -0.1, 0, lcout};
Point(21) = {lout, wout, 0, lcout};
Point(22) = {-lout, wout, 0, lcout};
Point(23) = {-lout, -0.1, 0, lcout};
Point(24) = {-lout, -0.6, 0, lcout};
Point(25) = {lout, -0.6, 0, lcout};
Line(20)={20,21};
Line(21)={21,22};
Line(22)={22,23};
Line(23)={23,20};
Line Loop(20) = {20, 21, 22, 23};
Physical Line("Outer") = {20, 21, 22, 23};
Line(23)={23,24};
Line(24)={24,25};
Line(25)={25,20};
Physical Line("Left") = {22,23};
Physical Line("Right") = {25,20};
Physical Line("Bottom") = {24};
Physical Line("Top") = {21};
Line Loop(20) = {20, 21, 22, 23, 24, 25};
/* Plane Surface */
Plane Surface(1) = {20, 10};
......
......@@ -6,7 +6,7 @@ import random
import csv
import matplotlib.pyplot as plt
filename = "tellurium_0_01/Drag.csv"
filename = "again81/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 = 1000
ranges = 100
print(data)
for i in range(len(react)):
if i <ranges:
......@@ -48,6 +48,6 @@ 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.ylim([-1,160])
plt.ylim([-1,16])
#plt.xlim([1.5,24])
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