3d testcase metzger et al. experiments

testcases/drop-3d/metzger/adapt/mesh.geo

+L = 0.05;
+H = 0.6;
+P = 0.02;
+y = 0;
+lc = 0.01;
+
+Point(1) = {-L, H, -P, lc};
+Point(2) = {-L, -H, -P, lc};
+Point(3) = {L, -H, -P, lc};
+Point(4) = {L, H, -P, lc};
+Point(5) = {-L, H, P, lc};
+Point(6) = {-L, -H, P, lc};
+Point(7) = {L, -H, P, lc};
+Point(8) = {L, H, P, lc};
+
+Line(1) = {1, 2};
+Line(2) = {2, 3};
+Line(3) = {3, 4};
+Line(4) = {4, 1};
+Line(5) = {5, 6};
+Line(6) = {6, 7};
+Line(7) = {7, 8};
+Line(8) = {8, 5};
+Line(9) = {5, 1};
+Line(10) = {4, 8};
+Line(11) = {2, 6};
+Line(12) = {3, 7};
+
+Line Loop(1) = {1:4};
+Line Loop(2) = {5:8};
+Line Loop(3) = {9,-4,10,8};
+Line Loop(4) = {2,12,-6,-11};
+Line Loop(5) = {12,7,-10,-3};
+Line Loop(6) = {1,11,-5,9};
+
+Plane Surface(1) = {1};
+Plane Surface(2) = {2};
+Plane Surface(3) = {3};
+Plane Surface(4) = {4};
+Plane Surface(5) = {5};
+Plane Surface(6) = {6};
+Surface Loop(1) = {1,-3,-2,4,1,-5,-2,6};
+Volume(1) = {1};
+Physical Surface("Z") = {1,2};
+Physical Surface("X") = {5,6};
+Physical Surface("Top") = {3};
+Physical Surface("Bottom") = {4};
+Physical Volume("Domain") = {1};
+Physical Point("PtFix") = {1};
+Mesh.Algorithm = 5;
+
+Merge "lc.pos";
+Field[1] = PostView;
+Field[1].IView = 0;
+Background Field = 1;
+
+Mesh.CharacteristicLengthExtendFromBoundary = 0;
+Mesh.CharacteristicLengthFromPoints = 0;

testcases/drop-3d/metzger/drop.py

+#!/usr/bin/env python
+from marblesbag import fluid3 as fluid
+from marblesbag import scontact3
+
+import numpy as np
+import os
+import time
+import shutil
+import random
+
+def genInitialPosition(filename, r, rout, rhop, compacity) :
+    p = scontact3.ParticleProblem()
+    p.load_msh_boundaries("mesh.msh", ["Top", "Bottom","X","Z"])
+    
+    N = compacity*rout**3/(r**3)
+    e = 2*rout/(6*N/np.pi)**(1./3.)
+    
+    for x in np.arange(rout, -rout, -e):
+        for y in np.arange(rout, -rout, -e):
+            for z in np.arange(rout, -rout, -e):
+                R2 = x**2 + y**2 + z**2
+                if R2<(rout-r)**2:
+                    p.add_particle((x+random.uniform(-e/2+r,e/2-r), y+random.uniform(-e/2+r,e/2-r), z+random.uniform(-e/2+r,e/2-r)), r, 4./3.* r**3 * np.pi * rhop);
+    
+    p.position()[:, 1] += 0.52
+    print(len(p.position()[:, 1] ))
+    p.write(filename)
+
+
+outputdir = "MetzgerB1"
+if not os.path.isdir(outputdir) :
+    os.makedirs(outputdir)
+
+t = 0
+ii = 0
+
+r=154e-6
+R = 3.3e-3
+compacity = 0.2
+drho = 35.4
+p = scontact3.ParticleProblem()
+#R = np.random.uniform(45e-06, 90e-06, len(x))
+
+#physical parameters
+g =  -9.81
+rhop = 2450
+nu = 1.17/1030.
+rho = 1030#rhop-drho/compacity
+V = 0.5 # todo : estimate V base on limit velocity
+print('V',V)
+tEnd = 300
+
+#numerical parameters
+lcmin = 0.0004 # approx r*100 but should match the mesh size
+dt = 5e-5
+alpha = 1e-3
+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")
+mu = nu*rho
+
+genInitialPosition(outputdir + "/part-00000", r, R, rhop, compacity)
+
+p = scontact3.ParticleProblem(outputdir+"/part-00000")
+
+print("r = %g, m = %g\n" %  (p.r()[0], p.mass()[0]))
+print("RHOP = %g" % rhop)
+outf = 1
+outf1 = 100000
+
+#strong_boundaries = [("Top",0,.0),("Top",1,0.),("Top",2,0.),("Box",1,0.),("Box",0,0.),("Box",2,0.),("Top",3,0.)]
+strong_boundaries = [("Top",0,.0),("Top",1,0.),("Top",2,0.),("Bottom",0,.0),("Bottom",1,0.),("Bottom",2,0.),("X",0,.0),("X",1,0.),("X",2,0.),("Z",0,.0),("Z",1,0.),("Z",2,0.),("Top",3,0.)]
+#strong_boundaries = [("Top",1,0.),("Bottom",1,0.),("X",0,.0),("Z",2,0.),("Top",3,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_vtk(outputdir,0,0)
+
+ii = 0
+jj = 0
+tEnd1 = 100
+#p.velocity()[:,1] = 2./9. *(compacity)* R**2*g*(rhop-rho)/mu
+print(p.velocity()[1,1])
+print(1-compacity)
+
+
+#for jj in range(3):
+#    dt1 = dt
+#    t = 0
+#    if jj!=0:
+#        fluid.adapt_mesh(0.01,50,1,50,4e-4,100000)
+#    while t<tEnd1:
+#        forces = fluid.compute_node_force(dt1)
+#        fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
+#        niter = fluid.implicit_euler(dt1)
+#        t += dt1
+#        if niter < 5 and dt1 < 20:
+#            dt1 *= 3
+#        print(dt1)
+
+#        if ii %outf == 0 :
+#            ioutput = int(ii/outf) + 1
+#            p.write_vtk(outputdir, ioutput, t)
+#            fluid.export_vtk(outputdir, t, ioutput)
+#        ii += 1
+
+ii = 0
+t = 0
+tic = time.clock()
+forces = g*p.mass()
+while t < tEnd : 
+    if (ii%5==0 and ii != 0):
+       fluid.adapt_mesh(0.01,50,1,50,4e-4,100000)
+       forces = fluid.compute_node_force(dt)
+    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/drop-3d/metzger/mesh.geo

+L = 0.05;
+H = 0.6;
+P = 0.02;
+y = 0;
+lc = 0.01;
+
+Point(1) = {-L, H, -P, lc};
+Point(2) = {-L, -H, -P, lc};
+Point(3) = {L, -H, -P, lc};
+Point(4) = {L, H, -P, lc};
+Point(5) = {-L, H, P, lc};
+Point(6) = {-L, -H, P, lc};
+Point(7) = {L, -H, P, lc};
+Point(8) = {L, H, P, lc};
+Point(9) = {0,H,0,lc};
+Point(10) = {0,-H,0,lc};
+
+Line(1) = {1, 2};
+Line(2) = {2, 3};
+Line(3) = {3, 4};
+Line(4) = {4, 1};
+Line(5) = {5, 6};
+Line(6) = {6, 7};
+Line(7) = {7, 8};
+Line(8) = {8, 5};
+Line(9) = {5, 1};
+Line(10) = {4, 8};
+Line(11) = {2, 6};
+Line(12) = {3, 7};
+
+Line Loop(1) = {1:4};
+Line Loop(2) = {5:8};
+Line Loop(3) = {9,-4,10,8};
+Line Loop(4) = {2,12,-6,-11};
+Line Loop(5) = {12,7,-10,-3};
+Line Loop(6) = {1,11,-5,9};
+
+Plane Surface(1) = {1};
+Plane Surface(2) = {2};
+Plane Surface(3) = {3};
+Plane Surface(4) = {4};
+Plane Surface(5) = {5};
+Plane Surface(6) = {6};
+Surface Loop(1) = {1,-3,-2,4,1,-5,-2,6};
+Volume(1) = {1};
+Physical Surface("Z") = {1,2};
+Physical Surface("X") = {5,6};
+Physical Surface("Top") = {3};
+Physical Surface("Bottom") = {4};
+Physical Volume("Domain") = {1};
+Physical Point("PtFix") = {1};
+
+
+Point(15) = {0,.52,0};
+Field[1] = Attractor;
+Field[1].NodesList = {15};
+
+Field[2] = Threshold;
+Field[2].DistMax = 0.02;
+Field[2].DistMin = 0.01;
+Field[2].LcMax = 0.01;
+Field[2].LcMin = 0.0008;
+Field[2].IField = 1;
+
+Background Field = 2;
+
+Mesh.CharacteristicLengthFromPoints = 0;
+Mesh.CharacteristicLengthExtendFromBoundary = 0;