Commit 76c05b81 authored by Matthieu Constant's avatar Matthieu Constant
Browse files

testcases polydisperse

parent 6b243e34
...@@ -74,7 +74,7 @@ forces = g*p.mass() ...@@ -74,7 +74,7 @@ forces = g*p.mass()
while t < tEnd : while t < tEnd :
# if (ii%50==0 and ii != 0): # if (ii%50==0 and ii != 0):
# fluid.adapt_mesh(0.01,50,4e-5,autoEpsilon) # fluid.adapt_mesh(0.01,50,4e-5,autoEpsilon)
forces = fluid.compute_node_force(dt,r*10) forces = fluid.compute_node_force(dt,r*100)
vn = p.velocity() + forces * dt / p.mass() vn = p.velocity() + forces * dt / p.mass()
vmax = np.max(np.hypot(vn[:, 0], vn[:, 1])) vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
nsub = max(1, int(np.ceil((vmax * dt * 4)/min(10*p.r())))) nsub = max(1, int(np.ceil((vmax * dt * 4)/min(10*p.r()))))
......
...@@ -8,18 +8,18 @@ import time ...@@ -8,18 +8,18 @@ import time
import shutil import shutil
import random import random
def genInitialPosition(filename, r, rout, rhop) : def genInitialPosition(filename, r, rhop) :
p = scontact2.ParticleProblem() p = scontact2.ParticleProblem()
p.load_msh_boundaries("mesh.msh", ["Top", "Bottom","Left","Right"]) p.load_msh_boundaries("mesh.msh", ["Top", "Bottom","Left","Right"])
p.add_particle((0, 0), r, r**2 * np.pi * rhop); p.add_particle((0, 0), r, r**2 * np.pi * rhop);
p.position()[:, 1] += 0.45 p.position()[:, 1] += 1.7
p.write(filename) p.write(filename)
outputdir = "output" outputdir = "outputc"
if not os.path.isdir(outputdir) : if not os.path.isdir(outputdir) :
os.makedirs(outputdir) os.makedirs(outputdir)
...@@ -51,7 +51,7 @@ shutil.copy("mesh.msh", outputdir +"/mesh.msh") ...@@ -51,7 +51,7 @@ shutil.copy("mesh.msh", outputdir +"/mesh.msh")
#scontact2Interface.MeshLoader(p, "funnel.msh", ("Funnel", "Top", "Bottom", "Lateral")) #scontact2Interface.MeshLoader(p, "funnel.msh", ("Funnel", "Top", "Bottom", "Lateral"))
p.write(outputdir+"/part-00000") p.write(outputdir+"/part-00000")
genInitialPosition(outputdir + "/part-00000", r, 1e-3, rhop) genInitialPosition(outputdir + "/part-00000", r, rhop)
p = scontact2.ParticleProblem(outputdir+"/part-00000") p = scontact2.ParticleProblem(outputdir+"/part-00000")
...@@ -72,7 +72,7 @@ ii = 0 ...@@ -72,7 +72,7 @@ ii = 0
tic = time.clock() tic = time.clock()
forces = g*p.mass() forces = g*p.mass()
while t < tEnd : while t < tEnd :
forces = fluid.compute_node_force(dt) forces = fluid.compute_node_force(dt,10*r)
vn = p.velocity() + forces * dt / p.mass() vn = p.velocity() + forces * dt / p.mass()
vmax = np.max(np.hypot(vn[:, 0], vn[:, 1])) vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
nsub = max(1, int(np.ceil((vmax * dt * 4)/min(10*p.r())))) nsub = max(1, int(np.ceil((vmax * dt * 4)/min(10*p.r()))))
......
L = 0.075; L = 0.075;
H = 2; H = 2;
y = 0; y = 0;
lc = 0.001; lc = 0.02;
Point(1) = {-L, H, 0}; Point(1) = {-L, H, 0, lc};
Point(2) = {-L, 0, 0}; Point(2) = {-L, 0, 0, lc};
Point(3) = {L, 0, 0}; Point(3) = {L, 0, 0, lc};
Point(4) = {L, H, 0}; Point(4) = {L, H, 0, lc};
Line(1) = {1, 2}; Line(1) = {1, 2};
Line(2) = {2, 3}; Line(2) = {2, 3};
Line(3) = {3, 4}; Line(3) = {3, 4};
Line(4) = {4, 1}; Line(4) = {4, 1};
Point(5) = {0,.25,0};
Point(6) = {0,H,0};
Line(5) = {5,6};
Line Loop(1) = {1:4}; Line Loop(1) = {1:4};
Plane Surface(1) = {1}; Plane Surface(1) = {1};
...@@ -24,19 +20,3 @@ Physical Line("Left") = {1}; ...@@ -24,19 +20,3 @@ Physical Line("Left") = {1};
Physical Line("Right") = {3}; Physical Line("Right") = {3};
Physical Surface("Domain") = {1}; Physical Surface("Domain") = {1};
Physical Point("PtFix") = {1}; Physical Point("PtFix") = {1};
Field[1] = Attractor;
Field[1].EdgesList = {5};
Field[1].NNodesByEdge = 200;
Field[2] = Threshold;
Field[2].DistMax = 0.08;
Field[2].DistMin = 0.04;
Field[2].LcMax = 0.08;
Field[2].LcMin = 0.05;
Field[2].IField = 1;
Background Field = 2;
Mesh.CharacteristicLengthFromPoints = 0;
Mesh.CharacteristicLengthExtendFromBoundary = 0;
...@@ -41,7 +41,7 @@ def genInitialPosition(filename, r1, r2, rout, rhop1, rhop2) : ...@@ -41,7 +41,7 @@ def genInitialPosition(filename, r1, r2, rout, rhop1, rhop2) :
p.write(filename) p.write(filename)
outputdir = "outputGoute3" outputdir = "outputGoutteStokesGoodNew"
if not os.path.isdir(outputdir) : if not os.path.isdir(outputdir) :
os.makedirs(outputdir) os.makedirs(outputdir)
...@@ -59,14 +59,14 @@ g = -9.81 ...@@ -59,14 +59,14 @@ g = -9.81
rho = 1200 rho = 1200
rhop1 = 2400 rhop1 = 2400
rhop2 = 2400 rhop2 = 2400
nu = 1e-4#/(10**0.5) nu = 1e-4#/(100**0.5)
V = 0.5 # todo : estimate V base on limit velocity V = 0.5 # todo : estimate V base on limit velocity
print('V',V) print('V',V)
tEnd = 1000 tEnd = 10
#numerical parameters #numerical parameters
lcmin = 0.00005 # approx r*100 but should match the mesh size lcmin = 0.00005 # approx r*100 but should match the mesh size
dt = 5e-3 dt = 5e-4
alpha = 2.5e-4 alpha = 2.5e-4
epsilon = alpha*lcmin**2 /nu#2e-2*c*h*2#2e-2*c*h # ?? not sure ??1e-5 epsilon = alpha*lcmin**2 /nu#2e-2*c*h*2#2e-2*c*h # ?? not sure ??1e-5
autoEpsilon = False autoEpsilon = False
...@@ -82,7 +82,7 @@ p = scontact2.ParticleProblem(outputdir+"/part-00000") ...@@ -82,7 +82,7 @@ p = scontact2.ParticleProblem(outputdir+"/part-00000")
print("r = %g, m = %g\n" % (p.r()[0], p.mass()[0])) print("r = %g, m = %g\n" % (p.r()[0], p.mass()[0]))
print("RHOP = %g" % rhop1) print("RHOP = %g" % rhop1)
outf = 40 outf = 50
outf1 = 100000 outf1 = 100000
strong_boundaries = [("Top",2,0.),("Top",1,0.),("Top",0,0.),("Bottom",1,0.),("Bottom",0,0.),("Left",0,0.),("Left",1,0.),("Right",0,0.),("Right",1,0.)] strong_boundaries = [("Top",2,0.),("Top",1,0.),("Top",0,0.),("Bottom",1,0.),("Bottom",0,0.),("Left",0,0.),("Left",1,0.),("Right",0,0.),("Right",1,0.)]
...@@ -98,9 +98,9 @@ ii = 0 ...@@ -98,9 +98,9 @@ ii = 0
tic = time.clock() tic = time.clock()
forces = g*p.mass() forces = g*p.mass()
while t < tEnd : while t < tEnd :
if (ii%100==0 and ii != 0) or ii==5: if (ii%100==0 and ii != 0) or ii==10:
fluid.adapt_mesh(0.01,100,1e-4,autoEpsilon) fluid.adapt_mesh(0.01,100,5e-5,autoEpsilon)
forces = fluid.compute_node_force(dt,0*r1) forces = fluid.compute_node_force(dt,10*r1)
vn = p.velocity() + forces * dt / p.mass() vn = p.velocity() + forces * dt / p.mass()
vmax = np.max(np.hypot(vn[:, 0], vn[:, 1])) vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
nsub = max(1, int(np.ceil((vmax * dt * 4)/min(10*p.r())))) nsub = max(1, int(np.ceil((vmax * dt * 4)/min(10*p.r()))))
......
...@@ -8,31 +8,32 @@ import time ...@@ -8,31 +8,32 @@ import time
import shutil import shutil
import random import random
def genInitialPosition(filename, r, rout, rhop) : def genInitialPosition(filename, r1, r2, rout, rhop1, rhop2) :
p = scontact2.ParticleProblem() p = scontact2.ParticleProblem()
p.load_msh_boundaries("mesh.msh", ["Top", "Bottom","Left","Right"]) p.load_msh_boundaries("mesh.msh", ["Top", "Bottom","Left","Right"])
x = np.arange(rout, -rout, -1.11e-4) e = 1.4e-3/sqrt(2)
y = np.arange(rout, -rout, -1.11e-4) x = np.arange(rout, -rout, -1.4e-3)
y = np.arange(rout, -rout, -1.4e-3)
x, y = np.meshgrid(x, y) x, y = np.meshgrid(x, y)
R2 = x**2 + y**2 R2 = x**2 + y**2
keep = R2 < (rout - r/2)**2 keep = R2 < (rout - r2/2)**2
x = x[keep] x = x[keep]
y = y[keep] y = y[keep]
N1 = 0 N1 = 0
N2 = 0 N2 = 0
for i in range(x.shape[0]) : for i in range(x.shape[0]) :
a = random.choice([1,2]) a = random.choice([1,2])
if a == 1 and N1 <= 128: if a == 1 and N1 <= 17:
N1 += 1 N1 += 1
p.add_particle((x[i]+random.uniform(-1.11e-4/3,1.11e-4/3), y[i]+random.uniform(-1.11e-4/3,1.11e-4/3)), r/2, r**2 * np.pi * rhop/4); p.add_particle((x[i]+random.uniform(-1.4e-3/3,1.4e-3/3), y[i]+random.uniform(-1.4e-3/3,1.4e-3/3)), r1, r1**2 * np.pi * rhop1);
elif a == 2 and N2 <= 128: elif a == 2 and N2 <= 17:
N2 += 1 N2 += 1
p.add_particle((x[i]+random.uniform(-1.11e-4/3,1.11e-4/3), y[i]+random.uniform(-1.11e-4/3,1.11e-4/3)), r, r**2 * np.pi * rhop); p.add_particle((x[i]+random.uniform(-1.4e-3/3,1.4e-3/3), y[i]+random.uniform(-1.4e-3/3,1.4e-3/3)), r2, r2**2 * np.pi * rhop2);
else: else:
if N1 > 128 and a == 1: if N1 > 17 and a == 1:
p.add_particle((x[i]+random.uniform(-1.11e-4/3,1.11e-4/3), y[i]+random.uniform(-1.11e-4/3,1.11e-4/3)), r, r**2 * np.pi * rhop); p.add_particle((x[i]+random.uniform(-1.4e-3/3,1.4e-3/3), y[i]+random.uniform(-1.4e-3/3,1.4e-3/3)), r2, r2**2 * np.pi * rhop2);
if N2 > 128 and a == 2: if N2 > 17 and a == 2:
p.add_particle((x[i]+random.uniform(-1.11e-4/3,1.11e-4/3), y[i]+random.uniform(-1.11e-4/3,1.11e-4/3)), r/2, r**2 * np.pi * rhop/4); p.add_particle((x[i]+random.uniform(-1.4e-3/3,1.4e-3/3), y[i]+random.uniform(-1.4e-3/3,1.4e-3/3)), r1, r1**2 * np.pi * rhop1);
print('N1',N1,'N2',N2,'N',N1+N2) print('N1',N1,'N2',N2,'N',N1+N2)
...@@ -41,21 +42,24 @@ def genInitialPosition(filename, r, rout, rhop) : ...@@ -41,21 +42,24 @@ def genInitialPosition(filename, r, rout, rhop) :
outputdir = "output" outputdir = "outputGoute"
if not os.path.isdir(outputdir) : if not os.path.isdir(outputdir) :
os.makedirs(outputdir) os.makedirs(outputdir)
t = 0 t = 0
ii = 0 ii = 0
r=0.36e-3/2 r1=0.36e-3/2/2
r2 = 0.78e-3/2/2
R = 4.1e-3/2
p = scontact2.ParticleProblem() p = scontact2.ParticleProblem()
#R = np.random.uniform(45e-06, 90e-06, len(x)) #R = np.random.uniform(45e-06, 90e-06, len(x))
#physical parameters #physical parameters
g = -9.81 g = -9.81
rho = 950 rho = 950
rhop = 2440 rhop1 = 2440
rhop2 = 2550
nu = 174e-6#/(10**0.5) nu = 174e-6#/(10**0.5)
V = 0.5 # todo : estimate V base on limit velocity V = 0.5 # todo : estimate V base on limit velocity
print('V',V) print('V',V)
...@@ -73,12 +77,12 @@ shutil.copy("mesh.msh", outputdir +"/mesh.msh") ...@@ -73,12 +77,12 @@ shutil.copy("mesh.msh", outputdir +"/mesh.msh")
#scontact2Interface.MeshLoader(p, "funnel.msh", ("Funnel", "Top", "Bottom", "Lateral")) #scontact2Interface.MeshLoader(p, "funnel.msh", ("Funnel", "Top", "Bottom", "Lateral"))
p.write(outputdir+"/part-00000") p.write(outputdir+"/part-00000")
genInitialPosition(outputdir + "/part-00000", r, 1e-3, rhop) genInitialPosition(outputdir + "/part-00000", r1, r2, R, rhop1, rhop2)
p = scontact2.ParticleProblem(outputdir+"/part-00000") p = scontact2.ParticleProblem(outputdir+"/part-00000")
print("r = %g, m = %g\n" % (p.r()[0], p.mass()[0])) print("r = %g, m = %g\n" % (p.r()[0], p.mass()[0]))
print("RHOP = %g" % rhop) print("RHOP = %g" % rhop1)
outf = 1 outf = 1
outf1 = 100000 outf1 = 100000
...@@ -96,7 +100,7 @@ forces = g*p.mass() ...@@ -96,7 +100,7 @@ forces = g*p.mass()
while t < tEnd : while t < tEnd :
# if (ii%50==0 and ii != 0): # if (ii%50==0 and ii != 0):
# fluid.adapt_mesh(0.01,50,4e-5,autoEpsilon) # fluid.adapt_mesh(0.01,50,4e-5,autoEpsilon)
forces = fluid.compute_node_force(dt) forces = fluid.compute_node_force(dt,1e-3)
vn = p.velocity() + forces * dt / p.mass() vn = p.velocity() + forces * dt / p.mass()
vmax = np.max(np.hypot(vn[:, 0], vn[:, 1])) vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
nsub = max(1, int(np.ceil((vmax * dt * 4)/min(10*p.r())))) nsub = max(1, int(np.ceil((vmax * dt * 4)/min(10*p.r()))))
......
import numpy as np
v0 = []
v1 = []
v2 = []
for k in range(0,1000):
with open("output/part-%05d" % k) as f :
for l in f :
w = l.split(" ")
if w[0] == "P" :
u = float(w[4])
print(u)
v0.append(u)
v0 = np.array(v0)
from matplotlib import pyplot as plt
plt.plot(v0)
plt.show()
Markdown is supported
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