dep2d+shaker

testcases/depot-2d/dep.py

+#!/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, r, ly, rhop) :
+    p = scontact2.ParticleProblem()
+    p.load_msh_boundaries("mesh.msh", ["Top", "Lateral","Bottom","BottomOut","TopOut"])
+    x = np.arange(-0.025+r, 0.025-r, r)
+    y = np.arange(0.06-r, .04, -r)
+    x, y = np.meshgrid(x, y)
+    x = x.flat
+    y = y.flat
+    for i in range(len(x)) :
+        rhop = random.choice([1100,1100,1100])
+        p.add_particle((x[i], y[i]), r/2, r**2/4 * np.pi * rhop);
+    p.write(filename)
+
+
+outputdir = "outputdep"
+if not os.path.isdir(outputdir) :
+    os.makedirs(outputdir)
+
+t = 0
+ii = 0
+
+r = 2e-4
+ly = 5e-2
+p = scontact2.ParticleProblem()
+#R = np.random.uniform(45e-06, 90e-06, len(x))
+
+#physical parameters
+g =  -9.81
+rho = 1000
+rhop = 1500
+nu = 1e-6
+V = 0.5 # todo : estimate V base on limit velocity
+print('V',V)
+tEnd = 100
+
+#numerical parameters
+lcmin = 0.001 # approx r*100 but should match the mesh size
+dt = 1e-4
+alpha = 2.5e-5
+epsilon = alpha*lcmin**2 /nu
+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", r, ly, rhop)
+
+p = scontact2.ParticleProblem(outputdir+"/part-00000")
+
+print("r = %g, m = %g\n" %  (p.r()[0], p.mass()[0]))
+print("RHOP = %g" % rhop)
+outf = 200
+outf1 = 100000
+
+ii = 0
+
+def outerBndV(x) :
+    bndLoop = 12000
+    if ii%bndLoop<11000 and ii%bndLoop>9500:
+        print(min(((ii-9500)%bndLoop)*5e-5,.05))
+        return min(((ii-9500)%bndLoop)*5e-5,.05)
+    else:
+        print(.05/((ii-11000)%bndLoop + 1))
+        return .05/((ii-11000)%bndLoop + 1)
+
+
+strong_boundaries = [("Top",2,0.),("TopOut",1,0),("Top",1,0),("BottomOut",1,0),("Bottom",1,0),("Lateral",0,0.),("Lateral",1,0.)]
+fluid = fluid.fluid_problem("mesh.msh",g,nu*rho,rho,epsilon,strong_boundaries)
+
+fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
+fluid.export(outputdir,0,0)
+
+
+tic = time.clock()
+forces = g*p.mass()
+while t < tEnd : 
+    forces = fluid.compute_node_force(dt)
+    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(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)
+        fluid.export_vtk(outputdir, t, ioutput)
+    ii += 1
+    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.clock() - tic))

testcases/depot-2d/mesh.geo

-L = 0.15;
-H = 2;
+L = 0.025;
+H = 0.06;
 y = 0;
-lc = 0.1;
+lc = 0.0015;
 
 Point(1) = {-L, H, 0, lc};
 Point(2) = {-L, -H, 0, lc};
 Point(3) = {L, -H, 0, lc};
 Point(4) = {L, H, 0, lc};
+
+Point(5) = {-8*L/9.,-H,0,lc};
+Point(6) = {-7*L/9.,-H,0,lc};
+Point(7) = {-6*L/9.,-H,0,lc};
+Point(8) = {-5*L/9.,-H,0,lc};
+Point(9) = {-4*L/9.,-H,0,lc};
+Point(10) = {-3*L/9.,-H,0,lc};
+Point(11) = {-2*L/9.,-H,0,lc};
+Point(12) = {-L/9.,-H,0,lc};
+
+Point(21) = {0.,-H,0,lc};
+Point(22) = {L/9.,-H,0,lc};
+Point(23) = {2*L/9.,-H,0,lc};
+Point(24) = {3*L/9.,-H,0,lc};
+Point(25) = {4*L/9.,-H,0,lc};
+Point(26) = {5*L/9.,-H,0,lc};
+Point(27) = {6*L/9.,-H,0,lc};
+Point(28) = {7*L/9.,-H,0,lc};
+Point(29) = {8*L/9.,-H,0,lc};
+
+Point(13) = {-8*L/9.,H,0,lc};
+Point(14) = {-7*L/9.,H,0,lc};
+Point(15) = {-6*L/9.,H,0,lc};
+Point(16) = {-5*L/9.,H,0,lc};
+Point(17) = {-4*L/9.,H,0,lc};
+Point(18) = {-3*L/9.,H,0,lc};
+Point(19) = {-2*L/9.,H,0,lc};
+Point(20) = {-L/9.,H,0,lc};
+
+Point(30) = {0.,H,0,lc};
+Point(31) = {L/9.,H,0,lc};
+Point(32) = {2*L/9.,H,0,lc};
+Point(33) = {3*L/9.,H,0,lc};
+Point(34) = {4*L/9.,H,0,lc};
+Point(35) = {5*L/9.,H,0,lc};
+Point(36) = {6*L/9.,H,0,lc};
+Point(37) = {7*L/9.,H,0,lc};
+Point(38) = {8*L/9.,H,0,lc};
+
 Line(1) = {1, 2};
-Line(2) = {2, 3};
-Line(3) = {3, 4};
-Line(4) = {4, 1};
-Line Loop(1) = {1:4};
+
+Line(2) = {2, 5};
+Line(3) = {5, 6};
+Line(4) = {6, 7};
+Line(5) = {7, 8};
+Line(6) = {8, 9};
+Line(7) = {9, 10};
+Line(8) = {10, 11};
+Line(9) = {11, 12};
+Line(10) = {12, 21};
+Line(11) = {21, 22};
+Line(12) = {22, 23};
+Line(13) = {23, 24};
+Line(14) = {24, 25};
+Line(15) = {25, 26};
+Line(16) = {26, 27};
+Line(17) = {27, 28};
+Line(18) = {28, 29};
+Line(19) = {29, 3};
+
+Line(20) = {3, 4};
+
+Line(21) = {4, 38};
+Line(22) = {38, 37};
+Line(23) = {37, 36};
+Line(24) = {36, 35};
+Line(25) = {35, 34};
+Line(26) = {34, 33};
+Line(27) = {33, 32};
+Line(28) = {32, 31};
+Line(29) = {31, 30};
+Line(30) = {30, 20};
+Line(31) = {20, 19};
+Line(32) = {19, 18};
+Line(33) = {18, 17};
+Line(34) = {17, 16};
+Line(35) = {16, 15};
+Line(36) = {15, 14};
+Line(37) = {14, 13};
+Line(38) = {13, 1};
+
+
+
+Line Loop(1) = {1:38};
 
 Plane Surface(1) = {1};
-Physical Line("Box") = {1,2,3};
-Physical Line("Top") = {4};
+Physical Line("Bottom") = {2,4,6,8,10,12,14,16,18};
+Physical Line("BottomOut") = {3,5,7,9,11,13,15,17,19};
+Physical Line("Lateral") = {1,20};
+Physical Line("Top") = {38,36,34,32,30,28,26,24,22};
+Physical Line("TopOut") = {21,23,25,27,29,31,33,35,37};
 Physical Surface("Domain") = {1};
 Physical Point("PtFix") = {1};

testcases/shaker/mesh.geo

 L = 0.025;
-H = 0.06;
+H = 0.025;
 y = 0;
 lc = 0.0015;