Commit 9b51f5e3 authored by Matthieu Constant's avatar Matthieu Constant
Browse files

testcase bidisperse

parent 96c74f23
L = 0.05;
L = 0.045;
H = 0.5;
y = 0;
lc = 0.1;
......
#!/usr/bin/env python
from marblesbag import fluid as fluid
from marblesbag import scontact2
import numpy as np
import os
import time
import shutil
import random
def genInitialPosition(filename, r1, r2, rout, rhop1, rhop2) :
p = scontact2.ParticleProblem()
p.load_msh_boundaries("mesh.msh", ["Top", "Bottom","Left","Right"])
e = .7e-4
x = np.arange(rout, -rout, -e)
y = np.arange(rout, -rout, -e)
x, y = np.meshgrid(x, y)
R2 = x**2 + y**2
keep = R2 < (rout - r2/2)**2
x = x[keep]
y = y[keep]
N1 = 0
N2 = 0
for i in range(x.shape[0]) :
a = random.choice([1,2])
if a == 1 and N1 <= 128:
N1 += 1
p.add_particle((x[i]+random.uniform(-e/3,e/3), y[i]+random.uniform(-e/3,e/3)), r1, r1**2 * np.pi * rhop1);
elif a == 2 and N2 <= 128:
N2 += 1
p.add_particle((x[i]+random.uniform(-e/3,e/3), y[i]+random.uniform(-e/3,e/3)), r2, r2**2 * np.pi * rhop2);
else:
if N1 > 128 and a == 1:
p.add_particle((x[i]+random.uniform(-e/3,e/3), y[i]+random.uniform(-e/3,e/3)), r2, r2**2 * np.pi * rhop2);
if N2 > 128 and a == 2:
p.add_particle((x[i]+random.uniform(-e/3,e/3), y[i]+random.uniform(-e/3,e/3)), r1, r1**2 * np.pi * rhop1);
print('N1',N1,'N2',N2,'N',N1+N2)
p.position()[:, 1] += 0.45
p.write(filename)
outputdir = "outputGoute3"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
t = 0
ii = 0
r1=25e-6/2
r2 = 25e-6/4
R = 1e-3
p = scontact2.ParticleProblem()
#R = np.random.uniform(45e-06, 90e-06, len(x))
#physical parameters
g = -9.81
rho = 1200
rhop1 = 2400
rhop2 = 2400
nu = 1e-4#/(10**0.5)
V = 0.5 # todo : estimate V base on limit velocity
print('V',V)
tEnd = 1000
#numerical parameters
lcmin = 0.00005 # approx r*100 but should match the mesh size
dt = 5e-3
alpha = 2.5e-4
epsilon = alpha*lcmin**2 /nu#2e-2*c*h*2#2e-2*c*h # ?? not sure ??1e-5
autoEpsilon = False
print('epsilon',epsilon)
shutil.copy("mesh.msh", outputdir +"/mesh.msh")
#scontact2Interface.MeshLoader(p, "funnel.msh", ("Funnel", "Top", "Bottom", "Lateral"))
p.write(outputdir+"/part-00000")
genInitialPosition(outputdir + "/part-00000", r1, r2, R, rhop1, rhop2)
p = scontact2.ParticleProblem(outputdir+"/part-00000")
print("r = %g, m = %g\n" % (p.r()[0], p.mass()[0]))
print("RHOP = %g" % rhop1)
outf = 40
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.)]
fluid = fluid.fluid_problem("mesh.msh",g,nu*rho,rho,epsilon,autoEpsilon,strong_boundaries)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.export(outputdir,0,0)
ii = 0
tic = time.clock()
forces = g*p.mass()
while t < tEnd :
if (ii%100==0 and ii != 0) or ii==5:
fluid.adapt_mesh(0.01,100,1e-4,autoEpsilon)
forces = fluid.compute_node_force(dt,0*r1)
vn = p.velocity() + forces * dt / p.mass()
vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
nsub = max(1, int(np.ceil((vmax * dt * 4)/min(10*p.r()))))
print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
for i in range(nsub) :
p.iterate(dt/nsub, forces)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.implicit_euler(dt)
t += dt
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
p.write_vtk(outputdir, ioutput, t)
p.write(outputdir+"/part-%05d" % ioutput)
p.write_boundary_vtk(outputdir, ioutput, t)
fluid.export(outputdir, t, ioutput)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.clock() - tic))
......@@ -33,7 +33,7 @@ Field[2] = Threshold;
Field[2].DistMax = 0.02;
Field[2].DistMin = 0.01;
Field[2].LcMax = 0.008;
Field[2].LcMin = 0.005;
Field[2].LcMin = 0.001;
Field[2].IField = 1;
Background Field = 2;
......
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